Solving Bessel Function using Runge Kutta

Question

I'm working on an assignment for a class of mine and I'm supposed to write a code using a program of my choice (I've chosen Matlab) to solve the Bessel function differential equation using the 4th order Runge-Kutta method. For reference the Bessel function DE is:

x^2*(J_n)''+x*(J_n)'+(x^2-n^2)*J_n=0.

I'm able to separate this into two coupled first order DEs by:

(J_n)'=Z_n and

(Z_n)'+(1/x)*Z_n+[(x^2-n^2)/x^2]*J_n=0.

I have no experience with Matlab nor any other programming language before this assignment. I know Matlab has the 'ode45' command but I'm supposed to write the code myself, not rely on Matlab's commands. So far I've been working on the n=0 case for the Bessel function but I keep getting an error when I try and plot the function. The current error I have says: "Undefined function or method 'J' for input arguments of type 'double'." But I don't know how to fix this error nor if my code is even correct. Could someone tell me where I've gone wrong or what is the correct way to write this code?

h=0.01;             %step size

J_0(1)=1;           %initial condition for J_0

Z_0(1)=1;           %initial condition for Z_0-This value should be zero
                    %but Matlab gives me an error. To fix this, I input
                    %Z_0(1)-1 to use the correct value for Z_0(1).

x(1)=0.001;         %first value of x

dZ(Z_0,J_0)=(-1/x)*(Z_0-1)-J_0;

for i=[1:1:10]

    dZ1=(-1/x)*(Z_0-1)-J_0;
    dJ1=(Z_0(1)-1)*h;

    dZ2=(-1/x)*(Z_0-1+0.5*h)-(J_0+0.5*h*dJ1);
    dJ2=((Z_0(1)-1)+dZ1)*h;

    dZ3=(-1/x)*(Z_0-1+0.5*h)-(J_0+0.5*h*dJ2);
    dJ3=((Z_0(1)-1)+dZ1+dZ2)*h;

    dZ4=(-1/x)*(Z_0-1+h)-(J_0+h*dJ3);
    dJ4=((Z_0(1)-1)+dZ1+dZ2+dZ3)*h;

    J(i+1)=J(i)+(h/6)*(dJ1+2*dJ2+2*dJ3+dJ4);
end

plot(J_0);

Thanks in advance for any help

Answer 1

Your problem is on the line:

J(i+1)=J(i)+(h/6)*(dJ1+2*dJ2+2*dJ3+dJ4);

In the right-hand side of your assignment operator you use the variable J that is never set before i is taking the value 1. Looks like a typo to me (should it be J_0 instead?)

Also, don't forget your index i when computing your dJ and dZ stuff in the for loop.