Problem with solving system of equations using symbolic toolbox

So my task is to figure out the constants named c1, c2 and c3 which are found in the f(r) function.

The relationship between these constants and my data is as follows:

We have the U1, U2, U3 and Ue functions.

We also have the L = U1+U2+U3-Ue;

U1, U2, U3, Ue and f(r) substituted into the L = U1+U2+U3-Ue function makes c1, c2 and c3 determinable by: @L/@c1 = 0 , @L/@c2 = 0, etc. (where @ makes for partial differentiation).

Now according to my problem and the way i am solving it, i should be getting a symbolic expression for c1,c2,c3 which i could use to suit the needs of my problem (that is my f(r) - a height function of a diaphragm- when set equal to say the given volume would produce the right results for the parameters I have an interest in).

The forum comes in to the part i am useless (and new) to this matlab environment. I give u the code this far:

syms r c1 c2 c3 V1 V2 z 

b = 0.008;
a = b/2;
hp = 0.0004;
hpzt = 0.0002;
hb = 0.00002;
%b,hp,hpzt,hb in (m)
vb = 0.3;
vp = 0.23;

e33 = 15.031*10^(-9);
%e33, e11 in (F/m)

s11 = 16.4 * 10^(-12);
s12 = -5.75 * 10^(-12);
%s11 and s12 in (m^2/N)
v = -s12/s11;
d31 = -171 * 10^(-12); 
%d31 in (C/N)
Ep = 160;
Eb = 5.17;
Epzt = 69;
%Ep,Eb in (GPa = 10^9 N/m^2)

neutralheight = 1/2 * (((Ep*hp^2)/(1-vp^2))+(Eb*(((hb+hp)^2)-hp^2))/(1-vb^2)+(Epzt*(((hpzt+hb+hp)^2)-((hp+hb)^2))/(1-v^2)))/(((Ep*hp^2)/(1-vp^2))+(Eb*hb/(1-vb^2))+(Epzt*hpzt/(1-v^2)));

h1 = hp+hb+(hpzt/2)-neutralheight;

f(r) = ((1-((4*r^2)/b^2))^2)*(c1+c2*(1-((4*r^2)/b^2)+c3*((1-((4*r^2)/b^2))^2)));

dfr = diff(f, r, 1);

df2r = diff(f, r, 2);

fi2 = (d31 / 4*d31^2 - 2*s11*e33*(1-v))*(df2r+dfr*(1/r));

fi1 = ((V1 - V2)/hp) - 2*h1*fi2;

U1 = ((pi*Ep*(((hp-neutralheight)^3)+ neutralheight^3))/3*(1-vp^2)) * int((r*df2r^2)+(2*vp*df2r*dfr)+(1/r)*dfr^2, r, 0, b/2);

U2 = ((pi*Eb*(((hp+hb-neutralheight)^3)+(hp-neutralheight)^3))/3*(1-vb^2))*int((r*df2r^2)+2*vb*df2r*dfr+(dfr^2)*(1/r), r, 0, a/2);

U3 = ((pi*(((h1+hpzt/2)^3)+(h1-hpzt/2)^3))/3*s11*(1-v^2)) * int((r*df2r^2)+2*v*df2r*dfr+(1/r)*dfr^2, r, 0, a/2) - ((2*pi*d31*(((h1+hpzt/2)^3)+(h1-hpzt/2)^3))/3*s11*(1-v)) * int((r*df2r + dfr) * fi2, r, 0, a/2) - (pi*d31*(((h1+hpzt/2)^2) + (h1-hpzt/2)^2)/2*s11*(1-v)) * int((r*df2r + dfr)*fi1, r, 0, a/2);

Ue = (pi*e33-((2*pi*d31^2)/(s11*(1-v))))* int(int(((fi1^2)+4*z*fi1*fi2 + 4*(z^2)*(fi2^2)*r), r, 0, a/2), z, h1-hpzt/2, h1+hpzt/2) + ((pi*d31)/(s11*(1-v))) * int(int((z*(fi1+fi2*2*z)*(df2r+(dfr/r))*r), r, 0, a/2), z, h1-hpzt/2, h1+hpzt/2);

L = U1+U2+U3-Ue;

f1 = diff(L, c1, 1);

f2 = diff(L, c2, 1);

f3 = diff(L, c3, 1);

This far, the result for f1 -which is c1 per se- comes as an expression that involves c2 and c3. The same goes for f2 and f3. How can i solve these three to get an answer for each one that doesnt involve the other two?

Thank you all in advance.