Double eigenvalue 
There are two cases:
- We can find two linearly independent eigenvectors

- We can find only one eigenvector
(deficient case). We need to find a "generalized eigenvector"
with
. This gives the solution
.
Case 1: there are two linearly independent eigenvectors
Example: consider the ODE system 
A = sym([2 0;0 2])
A =

Find the eigenvalues: solve the characteristic equation
This gives the eigenvalues
. r = eig(A) % find eigenvalues of A
r =

Find the eigenvectors: For each eigenvalue r find
with 
Here
and we can find two linearly independent eigenvectors, e.g.,
(actually, any two linearly independent vectors work). [V,D] = eig(A) % find matrix V with eigenvectors, diagonal matrix D with eigenvalues
V =

D =

Therefore we get a fundamental set of solutions
Solve the ODE system with dsolve:
syms y(t) [2 1] % define vector valued function y(t) = [y1(t);y2(t)]
sol = dsolve( diff(y)==A*y ); % find general solution
ysol = simplify( subs(y(t),sol) ) % in y(t) substitute the values from sol
ysol =

This is the general solution, with
. We obtain
by substituting
: Y1 = subs(ysol,{C1,C2},{1,0})
Y1 =

We obtain
by substituting
: Y2 = subs(ysol,{C1,C2},{0,1})
Y2 =

Plot vector field and trajectories in the phase plane
You need to download the file vectorfield.m and put it in the directory where your mlx files are. f = @(t,y) A*y; % define function f(t,y) for vectorfield
vectorfield(f,-20:2:20,-20:2:20); hold on
fplot(Y1(1),Y1(2),[-2,2],'b'); % plot trajectory for Y1 in blue
fplot(-Y1(1),-Y1(2),[-2,2],'b:'); % plot trajectory for -Y1 in blue, dotted
fplot(Y2(1),Y2(2),[-2,2],'k'); % plot trajectory for Y2 in black
fplot(-Y2(1),-Y2(2),[-2,2],'k:'); hold off; % plot trajectory for -Y2 in black, dotted
axis equal; axis([-20 20 -20 20]); xlabel('y_1'); ylabel('y_2')
title('four trajectories in the phase plane')
The general solution is
, so every trajectory is a straight line in the phase plane. Case 2: there is only one linearly independent eigenvector (DEFICIENT CASE)
Example: consider the ODE system 
A = sym([5 -4;1 1])
A =

Find the eigenvalues: solve the characteristic equation
This gives the eigenvalues
. r = eig(A) % find eigenvalues of A
r =

Find the eigenvectors: For each eigenvalue r find
with 
Here
and any solution of
has the form
with
. There is only one linearly independent eigenvector (deficient case). [V,D] = eig(A) % find matrix V with eigenvectors, diagonal matrix D with eigenvalues
V =

D =

Note that Matlab returns only one eigenvector in V.
This gives one solution
For a second linearly independent solution we need a "generalized eigenvector"
with
. We can find all eigenvectors and generalized eigenvectors with jordan(A) (find "Jordan normal form" of matrix A)
[V,D] = jordan(A) % find matrix V with eigenvectors and generalized eigenvectors
V =

D =

The matrix D is no longer diagonal. It contains "1" above the diagonal if the corresponding vector in V is not an eigenvector, but a generalized eigenvector. Here we get
- an eigenvector
satisfying
, yielding the solution 
- a generalized eigenvector
satisfying
, yielding the solution 
Solve the ODE system with dsolve:
syms y(t) [2 1] % define vector valued function y(t) = [y1(t);y2(t)]
sol = dsolve( diff(y)==A*y ); % find general solution
ysol = simplify( subs(y(t),sol) ) % in y(t) substitute the values from sol
ysol =

This is the general solution, with
. We obtain
by substituting
: Y1 = subs(ysol,{C1,C2},{1,0})
Y1 =

We obtain
by substituting
: Y2 = subs(ysol,{C1,C2},{0,1})
Y2 =

Plot vector field and trajectories in the phase plane
f = @(t,y) A*y; % define function f(t,y) for vectorfield
vectorfield(f,-1:.1:1,-1:.1:1); hold on
fplot(Y1(1),Y1(2),[-2,2],'b'); % plot trajectory for Y1 in blue
fplot(-Y1(1),-Y1(2),[-2,2],'b:'); % plot trajectory for -Y1 in blue, dotted
fplot(Y2(1),Y2(2),[-2,2],'k'); % plot trajectory for Y2 in black
fplot(-Y2(1),-Y2(2),[-2,2],'k:'); hold off; % plot trajectory for -Y2 in black, dotted
axis equal; axis([-1 1 -.5 .5]); xlabel('y_1'); ylabel('y_2')
title('four trajectories in the phase plane')