How to make my code run faster

This is the code I have written for solving a system of differential equations

%%CODE STARTS HERE

ti=0; %Initial time

tf=0.2; %Final time

%%Initial values of piston pressures(p1i-p9i)

p1i=0;

p2i=0;

p3i=0;

p4i=0;

p5i=0;

p6i=0;

p7i=0;

p8i=0;

p9i=0;

%%Initial value of discharge chamber pressure

p10i=0;

[t,P]=ode45('pistons_10',[ti tf],[p1i p2i p3i p4i p5i p6i p7i p8i p9i p10i]);

plot(t,P(:,10));

hold on;

%%CODE ENDS HERE

and the code for function "pistons_10" is

%%CODE STARTS HERE

function f = pistons_10(t,P)

d = 1.031*((10)^-2);  %Piston Diameter in meters

b = 0.335; %Angle of swash plate in radians

Rp = 0.021; %Piston pitch radius on barrel in meters

Vo = 1.534*((10)^-6); %Piston chamber initial volume in meter^3

w = 188.4; % Rotating speed in rad/s

B = 8.547*((10)^8); % Bulk modulus in pa

Cd1 = 0.675; % Discharge flow coefficient

Cd2 = 0.61; % Discharge flow coefficient

p = 860; % Fluid density in kg/m^3

phi = (2*pi)/9; % Phase Difference between pistons in radians

Av = 4.457*((10)^-6); % Discharge orifice area of valve

f=zeros(10,1);

D=zeros(10,1);

A=zeros(10,1);

E = (pi*((d)^2)*Rp*w*tan(b))/4;

C = (pi*((d)^2)*Rp*tan(b))/4;

z = (w*t)/pi;

for j=1:5

      Vc = 3.156*((10)^-(j-8)); %Discharge chamber control volume in meter^3

      for i=1:9

          if mod(floor(z+((2/9)*(i-1))),2)==0

              A(i)=pi*((d)^2)/4;

          else

              A(i)=0;

          end

          D(i)=Cd1*A(i)*sqrt(2/p);

          f(i,1)=(B/(Vo-(C*(1-cos((w*t)-((i-1)*phi)))))) * ((E*(abs(sin((w*t)-((i-1)*phi)))))-(D(i)*(sign(P(i)-P(10)))*(sqrt(abs(P(i)-P(10))))));

      end

          D(10) = Cd2*Av*(sqrt(2/p));

          f(10,1)=(B/Vc)* ((D(1)*(sign(P(1)-P(10)))*(sqrt(abs(P(1)-P(10))))) + (D(2)*(sign(P(2)-P(10)))*(sqrt(abs(P(2)-P(10))))) + (D(3)*(sign(P(3)-P(10)))*(sqrt(abs(P(3)-P(10))))) + (D(4)*(sign(P(4)-P(10)))*(sqrt(abs(P(4)-P(10))))) + (D(5)*(sign(P(5)-P(10)))*(sqrt(abs(P(5)-P(10))))) + (D(6)*(sign(P(6)-P(10)))*(sqrt(abs(P(6)-P(10))))) + (D(7)*(sign(P(7)-P(10)))*(sqrt(abs(P(7)-P(10))))) + (D(8)*(sign(P(8)-P(10)))*(sqrt(abs(P(8)-P(10))))) + (D(9)*(sign(P(9)-
  P(10)))*(sqrt(abs(P(9)-P(10))))) - (D(10)*(sqrt(abs(P(10))))));

end

end

%%CODE ENDS HERE

This code had been running past 24 hrs and matlab is still showing its status as 'busy'. Can someone suggest a way to make this code run faster

[Walter Roberson]

I have not determined why, but it appears that ode45 is quickly deciding that it needs a step size of about 1E-10 (or a small integer multiple of that), and then it ends up executing many many times.

It looks like f(7) is increasing at a rate of about 1200 per 1E-8 of time. Some of the others are a smaller slope but some of them are still more than half of that. That would reach roughly 2E10 by 0.2, I project, which is probably more than you are expecting. Your f(10) is the only one that is orders of magnitude smaller slope.

[Steven Lord]

Does the mathematical equation representing the ODE you're trying to solve really contain the sign function? [Please post the differential equations you're trying to solve -- not the code but the actual mathematical expressions.]