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
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
