I am a Matlab rookie. I need to solve numerically the following second order differential equations d^2x/dt^2 + w0_(el) * x = e/m_e * E - K3/m_e * x *y; d^2y/dt^2 + w0_(v) * y = - K_3/2M * x^2; I have started to deal with only the harmonic part of the problem. So I tried to solve d^2x/dt^2 + w0_(el) * x d^2y/dt^2 + w0_(v) * y with the following program t = 0:10^(-15):0.5*10^(-9); x0 = zeros(1,2); x0(1) = input('Insert the initial value of x'); x0(2) = input('Insert the initial value of dx/dt'); [t, x] = ode45(@harmonic, t, x0); plot(t, x(:, 1), 'g') title('Electronic position vs time'), xlabel('Time '), ylabel('Position') hold figure plot(t, x(:, 2), 'b') title('Nucleus position vs time'), xlabel('Time '), ylabel('Position') where the function "harmonic" is function ydot = harmonic(t, x) ydot = zeros(2,1); w_e = 10^30; w_v = 10^24; ydot(1) = x(3); ydot(2) = x(4); ydot(3) = -w_e*x(1); ydot(4) = -w_v*x(2); The outcome is a damping oscillation behaviour for x, which is of course meaningless because there aren't damping terms in the harmonic equations. So I am pretty desperate, if someone can help me I wuold be very glad. Thank you!