Function definitions are not permitted in this context.

I am trying to run a copied code, and I get that error. I named the file Sys2ODEsRK4

%%Defintion of the ODEs
dy1dt = @(t, y1, y2) y2;
dy2dt = @(t, y1, y2) (257600 - 5.1842*y1.^2) /(3000 - 80*t) - 32.2;
%%Solution Parameters
h = 1e-3; % Time Step
t_init = 0; % Initial Time
t_fin = 3; % Final Time
%%Set Initial Conditions (Rest)
y1 = 0;
dydt1 = 0;
%%Solve using RK4
[t, y, dydt] = Sys2ODEsRK4(dy1dt, dy2dt, t_init, t_fin, h, y1, dydt1);
%%Find Acceleration
d2ydt2 = (257600 - 5.1842*y.^2) ./ (3000 - 80*t) - 32.2;
%%Plot the solutions
subplot(3,1,1);
plot(t, y);
xlabel('Time (s)');
ylabel('Position (ft)');
subplot(3,1,2);
plot(t, dydt);
xlabel('Time (s)');
ylabel('Velocity (ft/s)');
subplot(3,1,3);
plot(t, d2ydt2);
xlabel('Time (s)');
ylabel('Accelerator (ft/s^2)');
%%The RK4 Method
function [t, x, y] = Sys2ODEsRK4(ODE1, ODE2, a, b, h, x1, y1)
% Sys2ODEsRK4 solves a system of two first-order initial value ODEs using
% fourth-order Runge-Kutta method.
% The independent variable is t, and the dependent variables are x and y.
% Input Variables:
% ODE1 Name for the function that calculates dx/dt
% ODE2 Name for the function that calculates dy/dt
% a The first value of t
% b The last value of t
% h The size of an increment
% x1 The initial value of x
% x2 The iniital value of yf
% Output Variables:
% t A vector with the t coordinate of the solution points.
% x A vector with the x coordinate of the solution points.
% y A vector with the y coordinate of the solution points.
t(1) = a; x(1) = x1; y(1) = y1;
n = (b-a)/h;
for i=1:n
t(i+1) = t(i) + h;
tm = t(i) + h/2;
Kx1 = ODE1(t(i), x(i), y(i));
Ky1 = ODE2(t(i), x(i), y(i));
Kx2 = ODE1(tm, x(i) + Kx1*h/2, y(i) + Ky1*h/2);
Ky2 = ODE2(tm, x(i) + Kx1*h/2, y(i) + Ky1*h/2);
Kx3 = ODE1(tm, x(i) + Kx2*h/2, y(i) + Ky2*h/2);
Ky3 = ODE2(tm, x(i) + Kx2*h/2, y(i) + Ky2*h/2);
Kx4 = ODE1(t(i+1), x(i) + Kx3*h, y(i) + Ky3*h);
Ky4 = ODE2(t(i+1), x(i) + Kx3*h, y(i) + Ky3*h);
x(i+1) = x(i) + (Kx1 + 2*Kx2 + 2*Kx3 + Kx4)*h/6;
y(i+1) = y(i) + (Ky1 + 2*Ky2 + 2*Ky3 + Ky4)*h/6;
end
end

Answer 1

[Walter Roberson]

It is only permitted to define functions inside a script from R2016b onwards.

In versions before that you would need to store the "function ... Sys2ODEsRK4" onward in file Sys2ODEsRK4.m; as you already used that name for the script, you would need to rename the script.

In R2016b and later you can define functions in the same file as scripts. However, none of the functions can have the same name as the script file has, and since your function is named Sys2ODEsRK4 you would need to rename the script to something other than Sys2ODEsRK4.m

[Mohamed afify]

The Version I am using is R2015S and I already named the script Sys2ODEsRK4.m

If I filled in the numbers through the command window it will work fine but it can't run with the variables inside the script

[Walter Roberson]

You have to store the below code in Sys2ODEsRK4.m

%%The RK4 Method
function [t, x, y] = Sys2ODEsRK4(ODE1, ODE2, a, b, h, x1, y1)
% Sys2ODEsRK4 solves a system of two first-order initial value ODEs using
% fourth-order Runge-Kutta method.
% The independent variable is t, and the dependent variables are x and y.
% Input Variables:
% ODE1 Name for the function that calculates dx/dt
% ODE2 Name for the function that calculates dy/dt
% a The first value of t
% b The last value of t
% h The size of an increment
% x1 The initial value of x
% x2 The initial value of yf
% Output Variables:
% t A vector with the t coordinate of the solution points.
% x A vector with the x coordinate of the solution points.
% y A vector with the y coordinate of the solution points.
t(1) = a; x(1) = x1; y(1) = y1;
n = (b-a)/h;
for i=1:n
t(i+1) = t(i) + h;
tm = t(i) + h/2;
Kx1 = ODE1(t(i), x(i), y(i));
Ky1 = ODE2(t(i), x(i), y(i));
Kx2 = ODE1(tm, x(i) + Kx1*h/2, y(i) + Ky1*h/2);
Ky2 = ODE2(tm, x(i) + Kx1*h/2, y(i) + Ky1*h/2);
Kx3 = ODE1(tm, x(i) + Kx2*h/2, y(i) + Ky2*h/2);
Ky3 = ODE2(tm, x(i) + Kx2*h/2, y(i) + Ky2*h/2);
Kx4 = ODE1(t(i+1), x(i) + Kx3*h, y(i) + Ky3*h);
Ky4 = ODE2(t(i+1), x(i) + Kx3*h, y(i) + Ky3*h);
x(i+1) = x(i) + (Kx1 + 2*Kx2 + 2*Kx3 + Kx4)*h/6;
y(i+1) = y(i) + (Ky1 + 2*Ky2 + 2*Ky3 + Ky4)*h/6;
end
end

You have to store the below in a file that is not named Sys2ODEsRK4.m . For example, you could name it Sys2ODEsRK4_driver.m

%%Definition of the ODEs
dy1dt = @(t, y1, y2) y2;
dy2dt = @(t, y1, y2) (257600 - 5.1842*y1.^2) /(3000 - 80*t) - 32.2;
%%Solution Parameters
h = 1e-3; % Time Step
t_init = 0; % Initial Time
t_fin = 3; % Final Time
%%Set Initial Conditions (Rest)
y1 = 0;
dydt1 = 0;
%%Solve using RK4
[t, y, dydt] = Sys2ODEsRK4(dy1dt, dy2dt, t_init, t_fin, h, y1, dydt1);
%%Find Acceleration
d2ydt2 = (257600 - 5.1842*y.^2) ./ (3000 - 80*t) - 32.2;
%%Plot the solutions
subplot(3,1,1);
plot(t, y);
xlabel('Time (s)');
ylabel('Position (ft)');
subplot(3,1,2);
plot(t, dydt);
xlabel('Time (s)');
ylabel('Velocity (ft/s)');
subplot(3,1,3);
plot(t, d2ydt2);
xlabel('Time (s)');
ylabel('Accelerator (ft/s^2)');