%define symbolic variables x, y
syms x y;
%define equations in two variables
circ = (x^2 + y^2 == 9); % a circle of radius 3
cub = (x^3 + 3*x == y);  % a cubic curve
%plot the two
clf;
ezplot(circ);
hold on;
ezplot(cub);
axis equal
%let's find the area below the cubic bounded by the circle.
%To do this, let's find the points of intersection
pts = solve(circ, cub)
%let's throw away complex solutions
%You can do it by hand (look at points with 0.000i)
%..or you can use indexing trickery, like I do here
real_x = sx((imag(pts.x)==0 & imag(sy)==0))
%(we really only need the x coordinates)
%Now we can set up an integral. We break the shape into two parts:
%below the cubic curve, and everything to the right

%The part below the cubic is bounded by these equations in x:
cub_x = x^3 + 3*x;
lower_circle = -sqrt(9-x^2);
%...and the bounds on x are given by intersections:
x1 = min(real_x);
x2 = max(real_x);
%so we get this part of the area:
part1 = int(cub_x - lower_circle, x, x1, x2);
part1 = double(part1) 
%The other part is everything to the right of x2
%This goes no further than 3, the radius of the circle
upper_circle = sqrt(9-x^2);
part2 = int(upper_circle - lower_circle, x, x2, 3);
part2 = double(part2)
%We get the area by adding two parts
area = part1 + part2
%However, all of the above could be avoided by noticing
%that the shape is symmetric around origin - so the area
%must be half of the area of the disk! We check this:
smarter_area = pi*(3^2)/2
%...and the answers are indeed the same.
