How to write code to solve Van der Pol oscillator using RK4 and Basic Euler?

I have been tasked to write a piece of code to solve; d2x/dt2−μ(1−x2)dx/dt+x=0 using the runge kutta 4th order and basic euler methods.

I have been having great difficulty with this so please will someone explain to me a) how my code (given further down) is unfit for purpose and how I could correct it or b) an alternative method to solving the problem.

In both cases please assume that I have very limited experience in Matlab...you won't be far wrong.

The problem statement is as follows; ____________________________________________________________

d2x/dt2−μ(1−x2)dx/dt+x=0

Use the implemented routines to find approximated solutions for the position of the oscillator in the interval 0 ≤ t ≤ 50 , with initial conditions x(0)=0.1, dx/dt(0)=0 when:

1. parameter μ = 0.1; solve by RK4 with Δt=0.1 and find a value of Δt for the Basic Euler, which guarantees that the solution does not diverge and the maximum difference between Basic Euler and RK4 over the last period is ≤ 10-2

Plot x(t) obtained by the Basic Euler method and the Runge-Kutta 4th order method and plot the difference of the solution given by the two methods over time.

I have attempted to codify a solution as detailed below; ____________________________________________________

%define ODEs

%x'=V

%V'=u(1-x^2)V-x

%define initial conditions

%let x,v describe the variables for the RK4 determination

%let X,V describe the variables for the Basic Euler determination.

x(1)=0.1;

X(1)=0.1;

v(1)=0;

V(1)=0;

u=0.1;

%define step size

h=0.1;

t=0:h:50;

F_xy

for j=1:10

     for i=1:(length(t)-1)

 t(i+1)=t(i)+h

 %Runge Kutta start

 %let k1 -> k4 be the k values for 'x'

 %let K1 -> K4 be the k values for 'v'

 k1=h*(v(i));

 K1=h*(u*(1-x(i)^2)*v(i)-x(i));

 k2=h*(v(i)+K1/2);

 K2=h*(u*(1-(x(i)+k1/2)^2)*(v(i)+K1/2)-(x(i)+k1/2));

 k3=h*(v(i)+K2/2);

 K3=h*(u*(1-(x(i)+k2/2)^2)*(v(i)+K2/2)-(x(i)+k2/2));

 k4=h*(v(i)+K3);

 K4=h*(u*(1-(x(i)+k3)^2)*(v(i)+K3)-(x(i)+k3));

 x(i+1)=x(i)+1/6*(k1+2*k2+2*k3+k4)

 v(i+1)=v(i)+1/6*(K1+2*K2+2*K3+K4)

 %Runge Kutta end

 %basic Euler start

 kX=h*(V(i));

 kV=h*(u*(1-X(i)^2)*V(i)-X(i));

 X(i+1)=X(i)+kX

 V(i+1)=V(i)+kV

 %Basic Euler end

 %define error function Ex

 Ex(i+1)=x(i+1)-X(i+1)

 error=abs(Ex)

     end

    if error(end)>0.01

        h=0.5*h

    else

        h=h

        break

    end

end

plot(t,x,t,Ex,'-o')

xlabel('time t')

ylabel('displacement x, error Ex')

figure;

___________________________________________________

When I run this it takes a stupendous amount of time to complete - I understand this is due to the 'for' loops calculating every tiny detail a million and one times but I don't know how to make a simpler, and more useful solution.

Thanks

[John BG]

Since this oscillator is quite common in 1st year electronics engineering modules, it's a popular case in MATLAB CENTRAL. Have a look:

http://www.mathworks.com/matlabcentral/fileexchange/?search_submit=fileexchange&query=van+der+pol&term=van+der+pol

John

[e-mail]