This is a MATLAB realization of your problem with matrix you suggest.
%% ;
k = 5; % number of nodes of the grid
D = sparse(1:k,1:k,2*ones(1,k),k,k); %diagonal
E = sparse(2:k,1:k-1,-1.25*ones(1,k-1),k,k); %bottom diagonal
A = (-1)*E + D + E'; %E' transposed bottom diagonal, i.e. upper one
x = linspace(0,4,k); % vector of independent variable
w_bound = zeros(1,k); w_bound(1) = 3/4; w_bound(k) = 1.25 * 24; % it's your w0 and wn I don't understant how you get it
b = -0.5 * (x'.^2 + 3) + w_bound'; % the right part of matrix equation
%% Constructing picture
figure(1)
hold on
plot(x, A \ b);
x_a = linspace(0,4,100);
y_a = (x_a + 1).^2 - exp(0.25 * x_a - 1) .* (1 / sin(sqrt(7))) .* sin(0.25 * sqrt(7) * x_a);
plot(x_a, y_a);
legend('your soulution', 'Amzoti''s Solution')
Picture with results of comparison your and Amzoti's solutions
It seems that your matrix A is not correct. Let's consider it again more accurately.
We have
$$y''=\frac{1}{2}y'-\frac{1}{2}y+\frac{x^2+3}{2}, ~~~~~y(0)=1, ~~y(4)=24,$$
and
$$y''=\frac{y_{i-1}-2y_i+y_{i+1}}{h^2}, y' = \frac{y_{i+1}-y_{i-1}}{2h}, y_0 = 1, y_n =24.$$
I used here central difference for the first derivative.
Let's sum up it together and collect terms:
$$ \frac{y_{i-1}-2y_i+y_{i+1}}{h^2} - \frac{1}{2}\frac{y_{i+1}-y_{i-1}}{2h} + \frac{1}{2}y_i = \frac{x_i^2+3}{2},y_0 = 1,y_n = 24$$, or
$$y_{i-1}\left(\frac{1}{h^2} + \frac{1}{4h}\right) + y_i\left(\frac{-2}{h^2}+\frac{1}{2}\right) + y_{i+1}\left(\frac{1}{h^2}-\frac{1}{4h}\right) = \frac{x_i^2+3}{2},y_0 = 1,y_n = 24$$.
I used the same $h = 1$ as you did.
MATLAB code:
%% the same steps only for my procedure
k = 5;
x = linspace(0, 4, k);
h = (x(end) - x(1)) / (k - 1);
D = sparse(1 : k, 1 : k,(-2 / h^2 + 0.5) * ones(1,k), k, k);
D_low = sparse(2 : k, 1 : k-1,(1 / h^2 + 1 / 4 /h) * ones(1, k - 1), k, k);
D_up = sparse(1 : k - 1, 2 : k,(1 / h^2 - 1 / 4 / h) * ones(1, k - 1), k, k);
A = D + D_low + D_up;
A(1,1) = 1; A(1,2) = 0;
A(end,end) = 1; A(end,end - 1) = 0;
b = (0.5 * (x.^2 + 3))'; b(1) = 1; b(end) = 24;
%% Pictures for finite difference solution
figure(2)
hold on
plot(x, A^(-1) * b);
x_a = linspace(0,4,100);
y_a = (x_a + 1).^2 - exp(0.25 * x_a - 1) .* (1 / sin(sqrt(7))) .* sin(0.25 * sqrt(7) * x_a);
plot(x_a, y_a);
legend('Finite difference method','Amzoti''s solution');
The matrix $A$ is:
$$\begin{bmatrix}
1 & 0 & 0 & & 0 & 0 \\
1.25 & -1.5 & 0.75 & \ddots & & 0 \\
0 & \ddots & \ddots & \ddots & \ddots & \\
& \ddots & \ddots & \ddots & \ddots & 0 \\
0 & & \ddots & 1.25 & -1.5 & 0.75 \\
0 & 0 & & 0 & 0 & 1
\end{bmatrix}.$$
Finite Difference and Amzoti's solutions picture
The condition for stability is that for eigenvalues $λ$ with negative real part you have $$|1+hλ|<1.$$ For your problem these eigenvalue are $λ=\frac12(-1\pm i\sqrt{3})$ and $λ=-\frac12(1+\sqrt5)$ so that the condition is
$$
(1-\tfrac12h)^2+\tfrac34h^2<1\iff h^2-h<0~~\text{ and }~~h(1+\sqrt5)<4
$$
and as $h>0$, the bound for stability is $h<\min(1, \frac{4}{1+\sqrt5})=1$.
However this only means that you avoid visibly, blatantly wrong behavior of the numerical solution, that is, that if the exact solution moves towards a fixed point, that the numerical solution were to explode away from it. The overall global error of size $O(h)$ or for variable time, $O((e^{Lt}-1)h)$, still remains. It may be advisable to take $h$ much smaller than the stability bound to get a sensible global error.
While the exact solution converges very rapidly towards the origin, with a bit more than one visible rotation about zero, the numerical solution with step size $h=0.5$ still somewhat shows that rapidity while $h=0.75$ has a rather slow convergence and $h=1$ looks undecided if it approaches the origin at all.
Best Answer
First off $\Sigma$ is of the wrong dimension. Since $C$ is a $3 \times 2$ matrix, so should $\Sigma$:
$$\Sigma=\left[\begin{matrix}1 & 0 \\ 0 & \sqrt 6 \\ 0 & 0\end{matrix}\right]$$
Next, you need to normalize $V$:
$$V=\left[\begin{matrix}-\frac{2}{\sqrt 5} & \frac{1}{\sqrt 5} \\ \frac{1}{\sqrt 5} & \frac{2}{\sqrt 5}\end{matrix}\right]$$
Now, we have:
$$U\Sigma=CV=\left[\begin{matrix}0 & \sqrt 5 \\ \frac{1}{\sqrt 5} & \frac{2}{\sqrt 5} \\ \frac{2}{\sqrt 5} & -\frac{1}{\sqrt 5}\end{matrix}\right]$$
Now, since the above is $U\Sigma$ and the first column of $\Sigma$ is $e_1$, the first column of this matrix must be the first column of $U$. Also, since the second column of $\Sigma$ is $\sqrt 6e_2$, the second column of $U$ must be the second column of $CV$ divided by $\sqrt 6$, so we get:
$$U=\left[\begin{matrix}0 & \sqrt{\frac 5 6} & ?? \\ \frac{1}{\sqrt 5} & \sqrt{\frac 2 {15}} & ?? \\ \frac{2}{\sqrt 5} & -\frac{1}{\sqrt 30} & ??\end{matrix}\right]$$
Finally, since $C$ has rank $2$ but $m=3$ rows, you need to pad the last column of $U$ with an orthogonal vector from the cokernel so that $U$ is an orthogonal matrix. Since the last column is orthogonal to the first two columns, it must be in the null space of the following matrix:
$$\left[\begin{matrix} 0 & \frac 1 {\sqrt 5} & \frac{2}{\sqrt 5} \\ \sqrt{\frac 5 6} & \sqrt{\frac 2 {15}} & -\frac{1}{\sqrt 30}\end{matrix}\right]$$
Basically, I just turned the two columns of $U$ into row vectors of the above matrix so that the null space of the above matrix is orthogonal to the columns of $U$. Now, find the RREF of the above matrix:
$$\left[\begin{matrix}1 & 0 & -1 \\ 0 & 1 & 2\end{matrix}\right]$$
From here, it should be clear that $(1,-2,1)^t$ spans the null space of this matrix. However, for $U$ to be orthogonal, we need this to be a unit vector, so normalize it by dividing by $\sqrt 6$:
$$\frac{1}{\sqrt{6}}\left[\begin{matrix}1 \\ -2 \\ 1\end{matrix}\right]=\left[\begin{matrix}\frac{1}{\sqrt 6} \\ -\sqrt{\frac 2 3} \\ \frac{1}{\sqrt 6}\end{matrix}\right]$$
Finally, we insert this column into $U$ to get:
$$U=\left[\begin{matrix}0 & \sqrt{\frac 5 6} & \frac{1}{\sqrt 6} \\ \frac{1}{\sqrt 5} & \sqrt{\frac 2 {15}} & -\sqrt{\frac 2 3} \\ \frac{2}{\sqrt 5} & -\frac{1}{\sqrt 30} & \frac{1}{\sqrt 6}\end{matrix}\right]$$