odefun=@(t,y)[-3*y(1)-y(2)+sin(t);y(1)-5*y(2)-2];
tspan=[0,10];
y0=[1,1];
h=0.01;
Nh=round((tspan(2)-tspan(1))/h);
[tn,un]=feuler(odefun,tspan,y0,Nh);

figure(1)
plot(tn,un(:,1),'r',tn,un(:,2),'b')
legend('y_1(t)','y_2(t)');
xlabel('t');
ylabel('numerical solution')
grid on