Forward Euler solution improvement

I would like to know if this code can be improved or is well-written as it is? it does not give errors. I'll appreciate any comments.Thanks

%part a

clear; clc;

% Parameters
alpha = 10;
beta_values = [2, 4];
h = 0.01;
t = 0:h:20;
N = length(t);

% Initial conditions
y1_0 = 0;
y2_0 = 2;

for b = 1:length(beta_values)
    beta = beta_values(b);

    
    y1 = zeros(1, N);
    y2 = zeros(1, N);

    % Initial values
    y1(1) = y1_0;
    y2(1) = y2_0;

    % Forward Euler method
    for n = 1:N-1
        f1 = alpha ...
             - y1(n) ...
             - (4*y1(n)*y2(n)) / (1 + y1(n)^2);

        f2 = beta * y1(n) * ...
             (1 - y2(n) / (1 + y1(n)^2));

        y1(n+1) = y1(n) + h * f1;
        y2(n+1) = y2(n) + h * f2;
    end

    figure;
    plot(t, y1, 'LineWidth', 1.5);
    xlabel('t');
    ylabel('y_1(t)');
    title(['Forward Euler: y_1 vs t, \beta = ', num2str(beta)]);
    grid on;

    figure;
    plot(y1, y2, 'LineWidth', 1.5);
    xlabel('y_1');
    ylabel('y_2');
    title(['Phase Portrait: y_2 vs y_1, \beta = ', num2str(beta)]);
    grid on;
end

%part b


% Parameters
alpha = 10;
beta_values = [3.4, 3.5, 3.6];   
h = 0.01;
t = 0:h:100;                    
N = length(t);

% Initial conditions
y1_0 = 0;
y2_0 = 2;

for b = 1:length(beta_values)
    beta = beta_values(b);

    
    y1 = zeros(1, N);
    y2 = zeros(1, N);

    % Initial values
    y1(1) = y1_0;
    y2(1) = y2_0;

    % Forward Euler method
    for n = 1:N-1
        f1 = alpha ...
             - y1(n) ...
             - (4*y1(n)*y2(n)) / (1 + y1(n)^2);

        f2 = beta * y1(n) * ...
             (1 - y2(n) / (1 + y1(n)^2));

        y1(n+1) = y1(n) + h * f1;
        y2(n+1) = y2(n) + h * f2;
    end

    
    figure;
    plot(t, y1, 'LineWidth', 1.5);
    xlabel('t');
    ylabel('y_1(t)');
    title(['y_1 vs t, \beta = ', num2str(beta)]);
    grid on;

    figure;
    plot(y1, y2, 'LineWidth', 1.5);
    xlabel('y_1');
    ylabel('y_2');
    title(['Phase Portrait y_2 vs y_1, \beta = ', num2str(beta)]);
    grid on;
end

Answer 1

[Alan Stevens]

It's clear and runs quickly. However, it's probably neater to pre-define your gradient functions like this:

% Gradient functions
f1 = @(alpha,y1,y2) alpha - y1 - (4*y1*y2) / (1 + y1^2);
f2 = @(beta,y1,y2)  beta * y1 * (1 - y2 / (1 + y1^2));

and then replace the forward Euler equations simply by:

    % Forward Euler method
    for n = 1:N-1
        y1(n+1) = y1(n) + h * f1(alpha,y1(n),y2(n));
        y2(n+1) = y2(n) + h * f2(beta,y1(n),y2(n));
    end

[Cesar]

Thanks, I would like to know if this code can also be improved? it works and does not give errors. Any comments will be appreciated.

clc; clear; close all;

% Parameters

T_final = 10;

u0 = [1; 0];

N_steps = [10000, 30000, 50000];

figure;

hold on;

colors = ['r', 'g', 'b'];

for i = 1:length(N_steps)

N = N_steps(i);

h = T_final / N;

t = 0:h:T_final;

u = zeros(2, length(t));

u(:, 1) = u0;

% Forward Euler

for n = 1:N

f = [u(2, n); -u(1, n)];

u(:, n+1) = u(:, n) + h * f;

end

plot(u(1, :), u(2, :), 'Color', colors(i), 'LineWidth', 1, ...

'DisplayName', ['N = ' num2str(N)]);

end

xlabel('u_1');

ylabel('u_2');

title('u_2 vs u_1 (Forward Euler)');

set(gcf, 'Position', [100, 100, 00, 500])

legend('show');

grid on;

axis equal;

[Cesar]

Also just wondering if oed45 is correctly implemented here?:

function orbit

mu = 0.012277471;

mu_hat = 1 - mu;

T_final = 17.06521656; % One full period

y0 = [0.994; 0; 0; -2.00158510637908252240537862224];

options = odeset('RelTol', 1e-12, 'AbsTol', 1e-12);

[t, y] = ode45(@(t, y) arenstorf_system(t, y, mu, mu_hat), [0 T_final], y0, options);

figure('Color', 'w', 'Position', [100, 100, 900, 700]);

plot(y(:,1), y(:,2), 'b--', 'LineWidth', 1.5);

hold on;

xlabel('u_1'); ylabel('u_2');

legend('ode45 solution');

grid on;

axis equal; % Ensures the orbit isn't visually distorted

end

function dy = arenstorf_system(~, y, mu, mu_hat)

u1 = y(1); u2 = y(2); v1 = y(3); v2 = y(4);

D1 = ((u1 + mu)^2 + u2^2)^(3/2);

D2 = ((u1 - mu_hat)^2 + u2^2)^(3/2);

du1 = v1;

du2 = v2;

dv1 = u1 + 2*v2 - mu_hat*(u1 + mu)/D1 - mu*(u1 - mu_hat)/D2;

dv2 = u2 - 2*v1 - mu_hat*u2/D1 - mu*u2/D2;

dy = [du1; du2; dv1; dv2];

end