I have a Hermite Cubic Finite Element Space on a computer in the form of Matlab m-files. More specifically, I can evaluate four "shape functions" $N_1, N_2, N_3,$ and $N_4$, for which the following holds on the boundary of the "reference element" $K_R= [-1,1]$: 
\begin{align*}
N_1(-1) = 1, N_1(1) = 0, N_1'(-1) = 0, N_1'(1) = 0,\\
N_2(-1) = 0, N_2(1) = 1, N_2'(-1) = 0, N_2'(1) = 0,\\
N_3(-1) = 0, N_3(1) = 0, N_3'(-1) = 1, N_3'(1) = 0,\\
N_4(-1) = 0, N_4(1) = 0, N_4'(-1) = 0, N_4'(1) = 1.
\end{align*}
Here is what they look like:
[![enter image description here][1]][1]


I can map an arbitrary element $K = [a,b]$ to $K_R$, evaluate the shape functions over $K_R$, then "map back", so that the above conditions on my shape functions hold on the boundary of arbitrary elements.

I want to use this Finite Element Space to approximate $u(x) = \sin(2\pi x)$ for $x \in \Omega = [2,3]$. To do so, I complete the following steps:

 1. Make an evenly spaced Mesh of the domain.
 2. Make relevant data structures needed to access the shape functions.
 3. Create an array of values of $u(x)$ evaluated at the Mesh nodes.
 4. Loop over each element: Extract the values of $u(x)$ at the Mesh node for this element $K$ with $u_1 = u(a)$ and $u_2 = u(b)$, where $a$ and $b$ are the endpoints of the element. Then do $u(x) \approx u_1 N_1(x) + u_2 N_2(x) + u_1 N_3(x) + u_2 N_4(x)$ for x = linspace(a,b,resolution). Then plot the result.

My problem is the Hermite Cubic Finite Element Space doesn't approximate $\sin(2 \pi x)$ very well using this method; the approximation wiggles a lot. Below is my plot using a Mesh with 100 elements, and using a resolution of 21. [![enter image description here][2]][2]

I repeated the same procedure as above except using a Lagrange Cubic Finite Element Space, so that the shape functions over $K_R$ satisfy
\begin{align*}
N_1(-1) = 1, N_1(-1/3) = 0, N_1(1/3) = 0, N_1(1) = 0, \\
N_2(-1) = 0, N_2(-1/3) = 1, N_2(1/3) = 0, N_2(1) = 0, \\
N_3(-1) = 0, N_3(-1/3) = 0, N_3(1/3) = 1, N_3(1) = 0, \\
N_4(-1) = 0, N_4(-1/3) = 0, N_4(1/3) = 0, N_4(1) = 1.
\end{align*}

Here is what the Lagrange shape functions look like:[![enter image description here][3]][3]


With $u_1, u_2, u_3, u_4$ corresponding to $-1, -1/3, 1/3, 1$, the approximation I used on each element was $u(x) \approx u_1 N_1(x) + u_2 N_2(x) +  u_3 N_3(x) +  u_4 N_4(x)$. Doing this gets me much nicer results:

[![enter image description here][4]][4]

What is the deal? Why does the Lagrange FE space work so nice, and the Hermite FE space fail? This might be more of a "numerical analysis" question rather than a "scientific computing" question, but I thought I would try here first.

____
EDIT: 

One of the comments suggested to use the formula
$u(x) \approx u_1 N_1(x) + u_2 N_2(x) + u_1' N_3(x) + u_2' N_4(x)$ instead of $u(x) \approx u_1 N_1(x) + u_2 N_2(x) + u_1 N_3(x) + u_2 N_4 (x)$. I agree that using this formula is what I should have done in this first place. I implement this by computing $u_1' = 2\pi \cos(2\pi a)$ and $u_2' = 2\pi \cos(2\pi b)$ for each element $K = [a,b]$. Then my Hermite approximation to $\sin(2 \pi x)$ looks like:

[![enter image description here][5]][5]

Which is clearly very wrong. In case it helps, here is the MATLAB code for my problem. I use a few of my own finite element computer programs here.

    % Make the mesh:
    domain = [2,3]; nElement = 100;
    Mesh = mesh_generator_1d(domain, nElement);
    
    % Form a data structure to access the finite element basis:
    iDegree = 3; Fem = femherm1d(Mesh, iDegree); hDegree = 1;
    
    % Create an array with the first 100 entries corresponding to the values
    % of sin(2*pi*x) evaluated at the Mesh nodes:
    uGlobal = sin(2*pi*Fem.point(1:nElement+1)); iDerivative = 0;
    
    % Extend the array so that the last 100 entries correspond to the values of
    % 2*pi*cos(2*pi*x) evaluated at the Mesh nodes:
    uGlobal(nElement+2:2*(nElement+1)) = 2*pi*cos(2*pi*Fem.point(nElement+2:2*(nElement+1)));
    
    % Set the resolution and prepare the figure:
    reso = 20; figure(1); clf; hold on;
    uInterpolate = zeros(1,(reso+1)*size(Mesh.element,2));
    
    for k = 1:size(Mesh.element,2) % Loop over each element in the mesh
        element = Mesh.node(Mesh.element(:,k)); % Extract element endpoints
        uLocal = uGlobal(Fem.T(:,k)); % uLocal(1) = u1, uLocal(2) = u2, uLocal(3) = u1', uLocal(4) = u2'
        x = linspace(element(1), element(2), reso+1); % x contains the points we are evaluating the shape functions at
        uInterpolate(((k-1)*(reso+1)+1):(k*(reso+1))) = evalfeherm1d(x, uLocal, ...
            element, hDegree, iDerivative); % Evaluate shape functions at x and scale by corresponding coefficients in uLocal 
        plot(x,uInterpolate(((k-1)*(reso+1)+1):(k*(reso+1)))); % Plot the result
    end
    hold off


  [1]: https://i.sstatic.net/u8BNv.png
  [2]: https://i.sstatic.net/DqMDi.png
  [3]: https://i.sstatic.net/ALWeh.png
  [4]: https://i.sstatic.net/J8GD9.png
  [5]: https://i.sstatic.net/3wr6l.png