RK4 help

Hi all,

Ive been trying to model one of J.C. Sprotts systems. The method being used is the 4th order runge-kutta method. The problem is I am unable to plot what I want to correctly. I know the plots Im getting are incorrect because the figures do not match up with the numbers (for example, output x1 at the bottom of the code).

I have tried using a hold on command and plotting inside my loop (t vs x1) and I have tried plotting outside of the loop. Neither works. Any ideas?

Thanks!

K.W.

Do not recommend prepackaged code please!

----------------

T = 0.1;
N = 100;

b = -5.6;
a = 0.5;

x1 = 0;
x2 = 0;
x3 = 0;

t = T*[0:N];

for n = 1:N

        u = sin(n*T);
        xdot11 = x2;
        xdot21 = x3;
        xdot31 = -x1 - b*(1-x1^2)*x2 - a*x3 + u;

        x1h = x1 + 0.5*T*xdot11;
        x2h = x2 + 0.5*T*xdot21;
        x3h = x3 + 0.5*T*xdot31;

        uh = sin((n+.5)*T);
        xdot12 = x2h;
        xdot22 = x3h;
        xdot32 = -x1h - b*(1-x1h^2)*x2h - a*x3h + uh;

        x1h = x1 + 0.5*T*xdot12;
        x2h = x2 + 0.5*T*xdot22;
        x3h = x3 + 0.5*T*xdot32;

        uh = sin((n+.5)*T);
        xdot13 = x2h;
        xdot23 = x3h;
        xdot33 = -x1h - b*(1-x1h^2)*x2h - a*x3h + uh;

        x1h = x1 + T*xdot13;
        x2h = x2 + T*xdot23;
        x3h = x3 + T*xdot33;

        uh = sin((n+1)*T);
        xdot14 = x2h;
        xdot24 = x3h;
        xdot34 = -x1h - b*(1-x1h^2)*x2h - a*x3h + uh;

        x1 = x1 + (T/6)*(xdot11 + 2*xdot12 + 2*xdot13 + xdot14)
        x2 = x2 + (T/6)*(xdot21 + 2*xdot22 + 2*xdot23 + xdot24);
        x3 = x3 + (T/6)*(xdot31 + 2*xdot32 + 2*xdot33 + xdot34);

end