Solve differential equation using Runge-Kutta

Hello everyone!

I have to solve the following equation by using the Runge-Kutta method:

all parameters are known, furthermore, Pin (function of t) and the timespan are vectors of dimensions 1x4096, so N(t) should be of the same dimension.

My code is the following:

n_dz=10;
dz=linspace(0,L,n_dz);

for i=2:n_dz
    Dz=dz(i)-dz(i-1)
    %%Compute N(t) using Runge-Kutte Method
    Fun=@(N) (I/(q*V))-A*N-B*N.^2-C*N.^3-Gamma*(S*(N-N_tran)...
        -k*(lambda-lambda_peak)^2).*(Pin/Dz)...
        *(exp(((Gamma*(S*(N-N_tran)-k*(lambda-lambda_peak)^2))-alfa_s)*dz(i))...
        -exp(((Gamma*(S*(N-N_tran)-k*(lambda-lambda_peak)^2))-alfa_s)*dz(i-1)))...
        /((Gamma*(S*(N-N_tran)-k*(lambda-lambda_peak)^2))-alfa_s)...
        *Dz/(V*h*f);
    y=Runge_Kutta_4(time,Dt,N_start,Fun);
    ...
    ...
end

with Runge_Kutta_4 implemented as:

function y=Runge_Kutta_4(timespan,h,N_start,fun)
    t=timespan;
    n=length(t);
    y=zeros(n,1);
    y(1)=N_start;
    for j=1:n-1
        k1 = h*feval(fun,t(j),y(j));
        k2 = h*feval(fun,t(j)+h/2,y(j)+k1/2);
        k3 = h*feval(fun,t(j)+h/2,y(j)+k2/2);
        k4 = h*feval(fun,t(j)+h,y(j)+k3);
        y(j+1) = y(j)+1/6*(k1+2*k2+2*k3+k4);
    end
end

Notice that in my code, the analysis over the total length "L" is divided into span of length "dz", therefore the last part of "Fun" is different from the formula shown above because it is due to an integral operation.

The following errors are shown:

    ??? Error using ==>
    @(N)(I/(q*V))-A*N-B*N.^2-C*N.^3-Gamma*(S*(N-N_tran)-k*(lambda-lambda_peak)^2).*(Pin/Dz)*(exp(((Gamma*(S*(N-N_tran)-k*(lambda-lambda_peak)^2))-alfa_s)*dz(i))-exp(((Gamma*(S*(N-N_tran)-k*(lambda-lambda_peak)^2))-alfa_s)*dz(i-1)))/((Gamma*(S*(N-N_tran)-k*(lambda-lambda_peak)^2))-alfa_s)*Dz/(V*h*f)
    Too many input arguments.

Error in ==> Runge_Kutta_4 at 7
        k1 = h*feval(fun,t(j),y(j));

Error in ==> SOA_v1 at 47
    y=Runge_Kutta_4(time,Dt,N_start,Fun);

Error in ==> RZ_source at 29

How can i solve this problem?

Thank you for your attention

Best regards

Answer 1

[Tim Berk]

I'm not sure I understand what your function (Fun) SHOULD be, but this is what you are telling it to be:

Fun = @(N) (I/(q*V))-A*N...

means Fun is a function of N, given by (I/(q*V))-A*N....

So you can put a value of N=9 into Fun using

feval(Fun,9)

And the output will be (I/(q*V))-A*9....

However, you are trying to put two values into the function by

k1 = h*feval(fun,t(j),y(j));

Hence the error Too many input arguments.