Problem solving a system of 4 ODE equations

beta = [1.875,4.694,7.855,10.996];

N1=[solve1(beta(1)),solve1(beta(2)),solve1(beta(3)),solve1(beta(4))];

N2=[solve2(beta(1),N1(1)),solve2(beta(2),N1(2)),solve2(beta(3),N1(3)),solve2(beta(4),N1(4))];

syms P Q l rho A E I x t real;

syms e(x,t) real;

syms W(x) real;

syms W1(x) W2(x) W3(x) W4(x) real;

W(x)=[W1(x),W2(x),W3(x),W4(x)];

syms omega(x,t) real;

rho=8030; % kg/m^3 %

A=4.81e-07; % m^2 %

E=200e09; % Pa %

I=7.75e-14; % N.sec/m^2 %

l=0.185; % m %

P=0.088; % N %

Q=0.244; % N %

W1(x)=solve3(N2(1),beta(1),N1(1)); W2(x)=solve3(N2(2),beta(2),N1(2)); W3(x)=solve3(N2(3),beta(3),N1(3)); W4(x)=solve3(N2(4),beta(4),N1(4));

syms phi1(t) phi2(t) phi3(t) phi4(t) real;

syms phi(t) real;

phi(t)=[phi1(t),phi2(t),phi3(t),phi4(t)];

omega(x,t)=sum(phi(t).*W(x));

syms d2phi1(t) d2phi2(t) d2phi3(t) d2phi4(t) real;

syms d2W1(x) d2W2(x) d2W3(x) d2W4(x) real;

syms d4W1(x) d4W2(x) d4W3(x) d4W4(x) real;

d2phi1(t)=diff(phi1(t),t,2); d2phi2(t)=diff(phi2(t),t,2); d2phi3(t)=diff(phi3(t),t,2); d2phi4(t)=diff(phi4(t),t,2);

d2W1(x)=diff(W1(x),x,2); d2W2(x)=diff(W2(x),x,2); d2W3(x)=diff(W3(x),x,2); d2W4(x)=diff(W4(x),x,2);

d4W1(x)=diff(W1(x),x,4); d4W2(x)=diff(W2(x),x,4); d4W3(x)=diff(W3(x),x,4); d4W4(x)=diff(W4(x),x,4);

e(x,t)=rho*A*((d2phi1(t)*W1(x))+(d2phi2(t)*W2(x))+(d2phi3(t)*W3(x))+(d2phi4(t)*W4(x)))+E*I*((phi1(t)*d4W1(x))+(phi2(t)*d4W2(x))+(phi3(t)*d4W3(x))+(phi4(t)*d4W4(x)))+P*((phi1(t)*d2W1(x))+(phi2(t)*d2W2(x))+(phi3(t)*d2W3(x))+(phi4(t)*d2W4(x)))-Q*dirac(x);

integrand=e(x,t).*W(x);

%Code for the first integral% %Splitting integral into parts%

expr1_part1=int(((W1(x))*rho*A*((d2phi1(t)*W1(x))+(d2phi2(t)*W2(x))+(d2phi3(t)*W3(x))+(d2phi4(t)*W4(x)))),x,0,l); expr1_part2=int((W1(x))*Q*dirac(x),x,0,l); expr1_part3=int(((W1(x))*P*((phi1(t)*d2W1(x))+(phi2(t)*d2W2(x))+(phi3(t)*d2W3(x))+(phi4(t)*d2W4(x)))),x,0,l); expr1_part4=int(((W1(x))*E*I*((phi1(t)*d4W1(x))+(phi2(t)*d4W2(x))+(phi3(t)*d4W3(x))+(phi4(t)*d4W4(x)))),x,0,l);

ode_eqn1=(vpa((expr1_part1-expr1_part2+expr1_part3+expr1_part4),4))==0;

%Code for the second integral% %Splitting integral into parts%

expr2_part1=vpa(int(((W2(x))*rho*A*((d2phi1(t)*W1(x))+(d2phi2(t)*W2(x))+(d2phi3(t)*W3(x))+(d2phi4(t)*W4(x)))),x,0,l),4); expr2_part2=vpa(int((W2(x))*Q*dirac(x),x,0,l),4); expr2_part3=vpa(int(((W2(x))*P*((phi1(t)*d2W1(x))+(phi2(t)*d2W2(x))+(phi3(t)*d2W3(x))+(phi4(t)*d2W4(x)))),x,0,l),4); expr2_part4=vpa(int(((W2(x))*E*I*((phi1(t)*d4W1(x))+(phi2(t)*d4W2(x))+(phi3(t)*d4W3(x))+(phi4(t)*d4W4(x)))),x,0,l),4);

ode_eqn2=(vpa((expr2_part1-expr2_part2+expr2_part3+expr2_part4),4))==0;

%Code for the third integral% %Splitting integral into parts%

expr3_part1=vpa(int(((W3(x))*rho*A*((d2phi1(t)*W1(x))+(d2phi2(t)*W2(x))+(d2phi3(t)*W3(x))+(d2phi4(t)*W4(x)))),x,0,l),4); expr3_part2=vpa(int((W3(x))*Q*dirac(x),x,0,l),4); expr3_part3=vpa(int(((W3(x))*P*((phi1(t)*d2W1(x))+(phi2(t)*d2W2(x))+(phi3(t)*d2W3(x))+(phi4(t)*d2W4(x)))),x,0,l),4); expr3_part4=vpa(int(((W3(x))*E*I*((phi1(t)*d4W1(x))+(phi2(t)*d4W2(x))+(phi3(t)*d4W3(x))+(phi4(t)*d4W4(x)))),x,0,l),4);

ode_eqn3=(vpa((expr3_part1-expr3_part2+expr3_part3+expr3_part4),4))==0;

%Code for the fourth integral% %Splitting integral into parts%

expr4_part1=vpa(int(((W4(x))*rho*A*((d2phi1(t)*W1(x))+(d2phi2(t)*W2(x))+(d2phi3(t)*W3(x))+(d2phi4(t)*W4(x)))),x,0,l),4); expr4_part2=vpa(int((W4(x))*Q*dirac(x),x,0,l),4); expr4_part3=vpa(int(((W4(x))*P*((phi1(t)*d2W1(x))+(phi2(t)*d2W2(x))+(phi3(t)*d2W3(x))+(phi4(t)*d2W4(x)))),x,0,l),4); expr4_part4=vpa(int(((W4(x))*E*I*((phi1(t)*d4W1(x))+(phi2(t)*d4W2(x))+(phi3(t)*d4W3(x))+(phi4(t)*d4W4(x)))),x,0,l),4);

ode_eqn4=(vpa((expr4_part1-expr4_part2+expr4_part3+expr4_part4),4))==0;

ode_eqns=[ode_eqn1;ode_eqn2;ode_eqn3;ode_eqn4];

S=dsolve(ode_eqns);

Dsolve returns complex functions for ph1(t), phi2(t), phi3(t) and phi4(t) when these functions are real valued. I would really appreciate a method to tell MATLAB that these functions are real or a good differential equation solver methhod in MATLAB. Essentially there are 4 ode equations that are obtained by solving 4 integrals and each equation is of the form c1*phi1(t)+c2*phi2(t)+c3*phi(t)+c4*phi4(t)+c5*d2phi1(t)+c6*d2phi2(t)+c7*d2phi3(t)+c8*d2phi4(t)=0 where d2phi#(t) is the second derivative of the phi functions.

[Elizabeth Reese]

You may try changing the IgnoreAnalyticConstraints to false in the call to dsolve. This changes the Algorithm as described at the link here.