function z = FEM1D(f, N)

%%%%%%%%%%%%%%%%%%
dx=0.1;
x=[0:dx:1]';
N=length(x);

%Aufbau der Matrizen
A=-2*diag(ones(N,1)); %Steifigkeitsmatrix
M=4*diag(ones(N,1));

%Nebendiagonalen befüllen
for i=1:N-1
    A(i,i+1)=1;
    A(i+1,i)=1;
    M(i,i+1)=1;
    M(i+1,i)=1;
end


%     a=0.5;
%     c=100;
%     f = c*(2*(2*a^2*c - 4*a*c*x + 2*c.*x.^2 - 1).*exp(-a^2*c + 2*a*c*x - c*x.^2));

% if f==@(x) (1-x)
%     f=(1-x);
% end

%f von runme
 %f=(1-x)

b=(M*f*dx^2/6);
% b=(f*dx^2/6);

%Randbedingungen
A(1,:)=0;
A(1,1)=1;
A(N,:)=0;
A(N,N)=1;

zL=0; %linker Rand
b(1)=zL;
zR=0; %rechter Rand
b(N)=zR;

%Lösung
z=A\b;

end