Is the ode45 function right?How to plot a graphic of 2 different function into one?

i got the fallowing code that is solving the superior order differential ecuation:

p2*d2_sigma+p1*d_sigma + p0*sigma=q2*d2_epsilon+q1*d_epsilon + q0*epsilon.

For a given epsilon(t) function, where t it's time,and q1, q2, q3, p1, p2, p3 are real parameters. I wanted to use ode45 function to find the sigma function using the initial conditions :

• sigma_initial=(q2/p2)*epsilon(0); % as in sigma_initial=sigma(0)
• dsigmadt_initial=(q1/p2-(q2*p1)/(p2^2))*epsilon(0)+(q2/p2)*d_epsilon(0);% as in dsigmadt_initial=d_sigma(0)

I want to ask if it's ok how i used ode45 to solve it. (i know it's something wrong because when i give the sigma function a constant it works right , but when i change it to a different function as an exponential or a trigonometric one, the sigma graphic does not goes to 0 after a give t interval as it should)

function serie_maxwell_kelvin

   clear all
   close all

   %ne definim functia epsilon  
  function [z] = epsilon(t)
  z=sin(t);
  end

   %definim derivatele functiei epsilon 
  syms t
  sym_epsilon = epsilon(t);
  clear d_epsilon d2_epsilon
  matlabFunction( diff(sym_epsilon, t), 'vars', t, 'file', 'd_epsilon.m');
  matlabFunction( diff(sym_epsilon, t, 2), 'vars', t, 'file', 'd2_epsilon.m');

   %introducem constantele(constanta de elasticitate si modulul lui Young)
  E1=input('E1=');
  niu1=input('niu1=');
  E2=input('E2=');
  niu2=input('niu2=');

   %calculam coeficienti pk si qk
  p2=niu2/E1;
  p1=1+E2/E1+niu2/niu1;
  p0=E2/niu1;
  q2=niu2;
  q1=E2;
  q0=0;

   %ne alegem intervalul de timp pe care vom reprezenta dependenta tensiuni
  t0=0;
  tf=100;
  N=101;
  t=linspace(t0,tf,101);

   %introducem conditiile de racordare a datelor initiale
  sigma_initial=(q2/p2)*epsilon(0);
  dsigmadt_initial=(q1/p2-(q2*p1)/(p2^2))*epsilon(0)+(q2/p2)*d_epsilon(0);
  %folosim functia ode45 ca sa ne rezolve ecuatia liniara de ordin superior
  %asociata sistemului de ecuatii diferentiale de ordinul intai
  [t,x]=ode45(@sistemul_cu_ec_de_gradul_1,t,[sigma_initial dsigmadt_initial]);

   %afisam graficul 
  figure(1)
  plot(t,x(:,1));
  xlabel('t');
  ylabel('sigma(t)');

   figure(2)
  plot(t,epsilon(t));
  %calculam sistemul de ecuatii diferentiale de ordinul intai asociat
  %ecuatiei liniare de ordin superior
  function[dxdt]=sistemul_cu_ec_de_gradul_1(t,x)
     dxdt_2=x(1);
     dxdt_1=-(p1/p2)*x(1)-(p0/p2)*x(2)+(q2/p2)*d2_epsilon(t)+(q1/p2)*d_epsilon(t)+(q0/p2)*epsilon(t);
     dxdt=[dxdt_1; dxdt_2];
  end

   end

Also i want to plot a graphic where i want as OX the epsilon values and on OY the sigma values on the same interval t. How can i do that?