Obtaining solution for non-homogeneous hyperbolic 1d PDE from homogeneous solution (Using explicit finite difference methods)

I'm having a problem with solving this non homogeneous hyperbolic 1d PDE:

I created a code that uses three explicit finite difference methods:

 if scheme == 1 % Lax Friedrichs   
    for n = 1:size(uu,2)-1
      for j = 2:size(uu,1)-1
            uu(j,n+1) = 0.5*(uu(j+1,n)+uu(j-1,n)) - 0.5*lambda*epsilon*(uu(j+1,n)-uu(j-1,n));
      end
    end
    
  elseif scheme == 2 % Lax Wendroff 
    for n = 1:size(uu,2)-1
      for j = 2:size(uu,1)-1
        uu(j,n+1) = uu(j,n) - 0.5*lambda*epsilon*(uu(j+1,n)-uu(j-1,n)) + 0.5*lambda^2*epsilon^2*(uu(j+1,n)-2*uu(j,n) + uu(j-1,n));
      end
      end 
        
  elseif scheme == 3 % Upwind
    
    for n = 1:size(uu,2)-1
      for j = 2:size(uu,1)-1
        uu(j,n+1) = uu(j,n) - 0.5*lambda*epsilon*(uu(j+1,n)-uu(j-1,n)) + 0.5*lambda*abs(epsilon)*(uu(j+1,n)-2*uu(j,n)+uu(j-1,n));
      end
    end
 end
 

But I think this code solves only the homogeneous part of the PDE, because when I plot the approximations and the function there is a clear error that it's probably due to the missing f(x,t) part of the approximation.

What should I do to adapt my code to make it work with a non homogeneous hyperbolic PDE?