Solving coupled second-order ODEs using ode45 (equations of motion)

Hello! I have browsed various threads on solving these types of problems and they've been very helpful. However, when I implemented the code, the answer makes no sense at all! I'm solving the equation where (Morse potential).

Here is my code:

syms M x(t) z(t) Y alpha V_0 l h_0

V = V_0 * ((exp(-alpha * z) - 1)^2 - 1)

% z motion
eqz = diff(z, 2) == -diff(V, z)/M
% x motion
eqx = diff(x, 2) == 0

[VF, Subs] = odeToVectorField(eqz, eqx)
ftotal = matlabFunction(VF, 'Vars', {t, Y, V_0, alpha, M})

% Parameters
A = 1E-10; % Angstrom
eV = 1.6E-19; % electron volt
amu = 1.66E-27;
% Depth of well
V_0 = 88 * eV/1000;
% Softness
alpha = 2/A;
% Mass
M = 3 * amu;
ps = 1E-12;

theta_i = 45;
phi_i = 0;

E_i = 50 * eV/1000; % Incident energy (J)
E_z = E_i * (cos(theta_i))^2;
p_mag = sqrt(2 * M * E_i); % Magnitude of incident momentum 
p_x_i = p_mag * sind(theta_i) * cosd(phi_i); % x-component of incident momentum
p_z_i = - p_mag * cosd(theta_i); % z-component of incident momentum

% initial conditions
ic = [20 * A; p_z_i/M; -10 * A; p_x_i/M];
tspan = [-2 * ps, 2 * ps];
[T,Y] = ode45(@(t,Y) ftotal(t, Y, V_0, alpha, M), tspan, ic)
% plot z against t
plot(T, Y(:, 1))

My problem is that the solution appears to just be a straight line, which doesn't make sense. Have I implemented ode45 in wrong?

Thanks for all the help!

Answer 1

[Wan Ji]

The odefun can be obtained using syms command

syms M x z Y alpha V_0 l h_0
V = V_0 * ((exp(-alpha * z) - 1)^2 - 1)
diff(V, z)/M

We get expressed by

-(2*V_0*alpha*exp(-alpha*z)*(exp(-alpha*z) - 1))/M

According to your explaining, the ode equations are

and

which can be implemented by an ode equation numerically

ftotal = @(t, Y, V_0, alpha, M)[ Y(2); % Y(1) = z; Y(2) = z';
    -(2*V_0*alpha*exp(-alpha*Y(1))*(exp(-alpha*Y(1)) - 1))/M;
    Y(4); % Y(3)=x; Y(4)=x';
    0;
    ];

Then you can use the your code following

ftotal = @(t, Y, V_0, alpha, M)[ Y(2); % Y(1) = z; Y(2) = z';
    -(2*V_0*alpha*exp(-alpha*Y(1))*(exp(-alpha*Y(1)) - 1))/M;
    Y(4); % Y(3)=x; Y(4)=x';
    0;
    ];
% Parameters
A = 1E-10; % Angstrom
eV = 1.6E-19; % electron volt
amu = 1.66E-27;
% Depth of well
V_0 = 88 * eV/1000;
% Softness
alpha = 2/A;
% Mass
M = 3 * amu;
ps = 1E-12;

theta_i = 45;
phi_i = 0;

E_i = 50 * eV/1000; % Incident energy (J)
E_z = E_i * (cos(theta_i))^2;
p_mag = sqrt(2 * M * E_i); % Magnitude of incident momentum 
p_x_i = p_mag * sind(theta_i) * cosd(phi_i); % x-component of incident momentum
p_z_i = - p_mag * cosd(theta_i); % z-component of incident momentum

% initial conditions
ic = [20 * A; p_z_i/M; -10 * A; p_x_i/M];
tspan = [-2 * ps, 2 * ps];
[T,Y] = ode45(@(t,Y) ftotal(t, Y, V_0, alpha, M), tspan, ic);
% plot z against t
plot(T, Y(:, 1))

The figure plotted is shown below

[Issac Jacob]

Thank you so much! This worked.

[Wan Ji]

@Issac Jacob I am happy to have my answer accepted by you! I need your support! ovo