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;