function [A, Mc, R, b, rv] = MyFirstPoissonAssembler(p);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% Syfte: Assemblerar matriser och vektorer för att, med cG(1)-metoden,
%% lösa differentialekvationen
%%
%% -(a(x)*u'(x))' + c(x)*u(x) = f(x),    x_min < x < x_max,
%%
%% med Robinrandvillkor
%%
%% a(x_min)*u'(x_min) = gamma(x_min)*(u(x_min) - g_D(x_min)) + g_N(x_min),
%% -a(x_max)*u'(x_max) = gamma(x_max)*(u(x_max) - g_D(x_max)) + g_N(x_max),
%%
%% Indata:
%% p:  mesh                  (Vektor med dimension 1 x N)
%%
%% Utdata:
%% A:  styvhetsmatris        (Matris med dimension N x N)
%% Mc: massmatris            (Matris med dimension N x N)
%% R:  randmatris            (Matris med dimension N x N)
%% b:  lastvektor            (Vektor med dimension N x 1)
%% rv: randvektor            (Vektor med dimension N x 1)

%% Anropade funktioner: a.m, c.m, f.m, robin_gamma.m, g_D.m, och g_N.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                     %
% Här börjar MyFirstPoissonAssembler. %
%                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N = length(p);  % Antal nodpunkter (inkl. randnoderna).

A = zeros(N, N);    % Initiera styvhetsmatris.
Mc = zeros(N, N);    % Initiera massmatris.
 R = zeros(N, N);    % Initiera randmatris.

b = zeros(N, 1);    % Initiera lastvektor.
rv = zeros(N, 1);    % Initiera randvektor.
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                      %
% Vi loopar nu över samtliga delintervall, och adderar deras bidrag    %
% till A, Mc, och b i tur och ordning. Detta kallas att assemblera.    %
%                                                                      %
% För att approximera integralerna över varje delintervall använder    %
% vi trapetsregeln för A och b, samt Simpsons formel för Mc.           %
%                                                                      %
% Notera att delintervall nummer i bidrar till följande matriselement  %
%                                                                      %
%    A(i, i),  A(i, i+1),  A(i+1, i),  A(i+1, i+1),                    %
%                                                                      %
%   Mc(i, i), Mc(i, i+1), Mc(i+1, i), Mc(i+1, i+1),                    %
%                                                                      %
% och följande vektorelement                                           %
%                                                                      %
%    b(i), b(i+1).                                                     %
%                                                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i = 1:(N-1)        % Loopa över delintervallen.

  xv = p(i); 	       % x-koordinat för delintervallets vänstra ändpunkt.
  xh = p(i+1);         % x-koordinat för delintervallets högra ändpunkt.
  xm = (xv + xh)/2;    % x-koordinat för delintervallets mittpunkt.
  
  h = xh - xv;         % Delintervallets längd.

  phi_v_xv = 1;        % Värdet i x = xv av den basfunktion, phi_v, som
                       % är lika med 1 i delintervallets vänstra ändpunkt.
  phi_v_xm = 1/2;      % Värdet i x = xm av phi_v.                     
  phi_v_xh = 0;        % Värdet i x = xh av phi_v.

  
  phi_h_xv = 0;        % Värdet i x = xv av den basfunktion, phi_h, som
                       %  är lika med 1 i delintervallets högra ändpunkt.
  phi_h_xm = 1/2;      % Värdet i x = xm av phi_h. 
  phi_h_xh = 1;        % Värdet i x = xh av phi_h.

  
  Dphi_v_xv = -1/h;    % Det konstanta värdet, d.v.s. speciellt värdet
                       % i x = xv, av derivatan (lutningen) av phi_v.
  
  Dphi_v_xh = -1/h;    % Värdet i x = xh av derivatan av phi_v.
  
  
  Dphi_h_xv = 1/h;     % Det konstanta värdet, d.v.s. speciellt värdet
                       % i x = xv, av derivatan av phi_h.
  
  Dphi_h_xh = 1/h;     % Värdet i x = xh av derivatan av phi_h.
  
  
  a_xv = a(xv);        % Värdet i x = xv av a(x) (definierad i filen a.m).
  a_xh = a(xh);        % Värdet i x = xh av a(x).
  
  c_xv = c(xv);        % Värdet i x = xv av c(x) (definierad i filen c.m).
  c_xm = c(xm);        % Värdet i x = xm av c(x).
  c_xh = c(xh);        % Värdet i x = xh av c(x).
  
  f_xv = f(xv);        % Värdet i x = xv av f(x) (definierad i filen f.m).
  f_xh = f(xh);        % Värdet i x = xh av f(x).
  
                       

  % Du får själv komplettera koden för assemblering av styvhetsmatrisen A!
  % Använd t.ex. trapetsregeln, eller integrera exakt.
  
  A(i, i) = A(i, i) + ??? % Integrand: a * Dphi_v * Dphi_v

  A(i, i+1) = A(i, i+1) + ??? % Integrand: a * Dphi_h * Dphi_v
  
  A(i+1, i) = A(i+1, i) + ??? % Integrand: a * Dphi_v * Dphi_h
  
  A(i+1, i+1) = A(i+1, i+1) + ??? % Integrand: a * Dphi_h * Dphi_h
  
  b(i)   = b(i)   + ??? % Integrand: f * phi_v
  b(i+1) = b(i+1) + ??? % Integrand: f * phi_h
  
  % Simpsons formel:
  % Kommentar: ... gör så att du kan fortsätta att skriva på nästa rad.

  Mc(i, i)     = Mc(i, i)     + ... % Integrand: c * phi_v * phi_v
      (c_xv*phi_v_xv*phi_v_xv + 4*c_xm*phi_v_xm*phi_v_xm + c_xh*phi_v_xh*phi_v_xh)*h/6;
  Mc(i, i+1)   = Mc(i, i+1)   + ... % Integrand: c * phi_h * phi_v
      (c_xv*phi_h_xv*phi_v_xv + 4*c_xm*phi_h_xm*phi_v_xm + c_xh*phi_h_xh*phi_v_xh)*h/6;
  Mc(i+1, i)   = Mc(i+1, i)   + ... % Integrand: c * phi_v * phi_h
      (c_xv*phi_v_xv*phi_h_xv + 4*c_xm*phi_v_xm*phi_h_xm + c_xh*phi_v_xh*phi_h_xh)*h/6;
  Mc(i+1, i+1) = Mc(i+1, i+1) + ... % Integrand: c * phi_h * phi_h
      (c_xv*phi_h_xv*phi_h_xv + 4*c_xm*phi_h_xm*phi_h_xm + c_xh*phi_h_xh*phi_h_xh)*h/6;
  
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                     %
% Vi adderar nu randbidragen från x = x_min och x = x_max.            %
%                                                                     %  
% Nodpunkten x = x_min bidrar till följande element i matrisen 'R'    %
%                                                                     %
%   R(1, 1),                                                          %
%                                                                     %
% och till följande element i vektorn 'rv'                            %
%                                                                     %
%   rv(1).                                                            % 
%                                                                     %
% Nodpunkten x = x_max bidrar till följande element i 'R'             %
%                                                                     %
%   R(N, N),                                                          %
%                                                                     %
% och till följande element i 'rv'                                    %
%                                                                     %
%   rv(N).                                                            % 
%                                                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Bidrag från nodpunkten x = x_min:

R(1, 1) = robin_gamma(1);    % x_min har id-nummer 1

rv(1)   = ???   % Du får själv komplettera koden för beräkning av rv(1)!


% Bidrag från nodpunkten x = x_max:

R(N, N) = ???   % x_max har id-nummer 2

rv(N)   = ???   % Du får själv komplettera koden för beräkning
                % av R(N,N) och rv(N)!


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                  %
% Slut på MyFirstPoissonAssembler. %
%                                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
