tmin = -1; tmax = 5;
ymin = -3; ymax = 8;
dt = 0.2; dy = 0.2;
 
[t, y] = meshgrid(tmin:dt:tmax, ymin:dy:ymax);
 
yp = 2*sin(3*t) - y;
 
% Normalize
L = sqrt(1 + yp.^2);
quiver(t, y, 1./L, yp./L, 0.5, 'k');  % direction field (black)
hold on;
 
% y(t) = (1/5)*(sin(3t) - 3*cos(3t)) + C*exp(t)
y_general = @(t, C) (1/5)*(sin(3*t) - 3*cos(3*t)) + C.*exp(-t);
 
C0 = (3/5); % y(0)=0
C1 = (8/5); % y(0)=1
C2 = (13/5); % y(0)=2
 
fplot(@(t) y_general(t, C0), [tmin tmax], 2);
fplot(@(t) y_general(t, C1), [tmin tmax], 2);
fplot(@(t) y_general(t, C2), [tmin tmax], 2);
 
axis([tmin tmax ymin ymax]);
grid on;
hold off;