function u=solve(coordinates,elements,dirichlet,a,b,c,f,u_d)

[A,B,C,rhs] = assemble(coordinates,elements,a,b,c,f);

M=A+B+C;

% Dirichlet Bedingungen
u = zeros(size(A,2),1);
u(dirichlet)=u_d;

% Modifiziere rechte Seite für Dirichlet Daten
rhs=rhs-M*u;

% Loese das LGS nur an den inneren Knoten 
freenodes=setdiff(1:size(elements,1)+1,dirichlet); 
u(freenodes)=M(freenodes,freenodes)\rhs(freenodes);