function [dx,ff1,ff2]=myfun(t,x)
t
beita=26;
mn=0.004;
z1=46;
z2=43;
z3=122;
T_in=200;
T_out=80;
roug1=7.8E3;
roug2=7.8E3;
roug3=7.8E3;
alphan=20;
alphat=atand(tand(alphan)/cosd(beita));
d1=z1*mn/cosd(beita)/1000;
db1=d1*cosd(alphat)/1000;
d2=z2*mn/cosd(beita)/1000;
db2=d2*cosd(alphat)/1000;
d3=z3*mn/cosd(beita)/1000;
db3=d3*cosd(alphat)/1000;
bp1=116/1000;
bp2=116/1000;
bp3=116/1000;
bp=116/1000;
I1=((roug1*pi*(d1/2)^2*bp1)*(d1/2)^2)/2;
I2=((roug2*pi*(d2/2)^2*bp2)*(d2/2)^2)/2;
I3=((roug3*pi*(d3/2)^2*bp3)*(d1/2)^2)/2;
m1=roug1*pi*(d1/2)^2*bp1;
m2=roug2*pi*(d2/2)^2*bp2;
m3=roug3*pi*