clear all clc syms y(t) dy = diff(y) deq = diff(y, 2) + dy + y == 0 soln = dsolve(deq, y(0) == 1, dy(0) == 0) y_eval = double(vpa(subs(soln, t, 4)))