second order variable coefficients ODE -so confused !!!

Hi matlab person:

    Nowdays, I am solving a second order variable coefficients ODE {x^2 y''+x y'-sin(y)cos(y)+x sin(y)^2-x^2 sin(y)cos(y)=0}. I feel there is no analytical solution, so I want to use ODE to get numerical solution.But I don't know how I should solve second order variable coefficients ODE through ODE45.

     Any suggestion, any help is ok.

    Thanks very much.

equation: x^2 y''+x y'-sin(y)cos(y)+x sin(y)^2-x^2 sin(y)cos(y)=0

boundary condition: y(0)=0;y'(0)=-0.126

[Yash]

CHECK MY RUNGE KUTTA CODE FOR THAT, MIGHT HELP, ITS FOR TWO VARIABLES

Following matlab functions will be used. function dxdt = f(t,x,y) % main function dxdt =2*y+x+t;

function dydt = g(t,x,y) % main function dydt =sin(x)+exp(t);

function [tout,yout,xout] = rk2variable(f,g,tspan,y0,N,x0) h = (tspan(2)-tspan(1))/N; t = tspan(1); tout = t; y = y0(:); yout = y.'; x = x0(:); xout = x.'; for n=1:N k1 =feval(f,t,x,y); j1=feval(g,t,x,y); k2=feval(f,t+h/2, x + (h/2)*k1, y + (h/2)*j1); j2=feval(g, t + h/2, x + (h/2)*k1, y + (h/2)*j1); k3=feval(f, t + h/2, x + (h/2)*k2, y + (h/2)*j2); j3=feval(g, t + h/2, x + (h/2)*k2, y + (h/2)*j2); k4=feval(f, t + h, x + h* k3, y + h* j3); j4=feval(g, t + h, x + h* k3, y + h* j3); x = x + (h/6)*(k1 + 2*k2 + 2*k3 + k4); y = y + (h/6)*(j1 + 2*j2 + 2*j3 + j4); t = t+h; yout = [yout; y.']; tout = [tout; t]; xout = [xout; x.']; end

This m file will be used to compute the values at different values of h.

clear tmin = 0; tmax = 1; y0 = 0; h =0.01; N = round((tmax-tmin)/h); x0=0; [tout,yout,xout] = rk2variable(@f,@g,[tmin,tmax],y0,N,x0); plot(tout,yout,tout,xout) legend('Y','X')

[Yash]

Following matlab functions will be used.
function dxdt = f(t,x,y) % main function
dxdt =2*y+x+t;

function dydt = g(t,x,y) % main function
dydt =sin(x)+exp(t);

function [tout,yout,xout] = rk2variable(f,g,tspan,y0,N,x0)
h = (tspan(2)-tspan(1))/N;
t = tspan(1); tout = t;
y = y0(:); yout = y.';
x = x0(:); xout = x.';
for n=1:N
k1 =feval(f,t,x,y);
j1=feval(g,t,x,y);
k2=feval(f,t+h/2, x + (h/2)*k1, y + (h/2)*j1);
j2=feval(g, t + h/2, x + (h/2)*k1, y + (h/2)*j1);
k3=feval(f, t + h/2, x + (h/2)*k2, y + (h/2)*j2);
j3=feval(g, t + h/2, x + (h/2)*k2, y + (h/2)*j2);
k4=feval(f, t + h, x + h* k3, y + h* j3);
j4=feval(g, t + h, x + h* k3, y + h* j3);
x = x + (h/6)*(k1 + 2*k2 + 2*k3 + k4);
y = y + (h/6)*(j1 + 2*j2 + 2*j3 + j4);
t = t+h;
yout = [yout; y.']; tout = [tout; t];
xout = [xout; x.'];
end

This m file will be used to compute the values at different values of h.

clear
tmin = 0; tmax = 1; y0 = 0;  h =0.01;
N = round((tmax-tmin)/h);
x0=0;
[tout,yout,xout] = rk2variable(@f,@g,[tmin,tmax],y0,N,x0);
plot(tout,yout,tout,xout)
legend('Y','X')

[Jan]

@Yash: I definitely would not hand-code a Runge-Kutta manually, when Matlab offers more powerful ODE functions.

[Jan]

@petter: Please consider, that the term "urgent" is not apropriate, when the answers are craeted by volunteers. Even the "please" does not shift the tone back to a friendly level sufficently.

Answer 1

[Jan]

You can convert the ODE system of order k for a problem of dimension d into a system of order 1 and dimension k*d. Instructions can be found at Wikipedia or in your Numerics script.

Then the documentation of ODE45 will help you to implement this in Matlab, especially the examples on "doc ode45".

[petter]

Jan simon, thank you very much for your reply and for your suggestion. I have thought this problem several days and asked other persons. This problem comes from the paper which I am preparing.So I used the term"urgent",may be that is not appropriate. sorry for this.

The following is my code, but there are several problem:

first defining function ***********************************************

function yprime = skyrmion(x,y)

k=0.95;

yprime = [y(2); -1/x*y(2)+1/x^2*sin(y(1))*cos(y(1))-4*k/pi* (1/x)*sin(y(1))^2+sin(y(1))*cos(y(1))]; ***********************************************************

Then transfer function

*******************************************************

xspan=[0.01,10];

y0=[-0.126;pi];

[x,y]=ode45('skyrmion',xspan,y0);

plot(x(:),y(:,2))

******************************************************

problem one:

Because The ODE was divided by x^2, so the x span can't use x==0. whether the code can be modified to not divide x^2.

problem two:

the results get by using my code are not consist with some papers.

so I want to know my code is right or wrong.

Thank you again.