Beam Equation Using Singularity Functions
Beam loading can occur in the form of singularity functions. By singularity functions we mean functions with discontinuity in either the slope, or the value of the function at one or more places. As discussed by Prof. S. Banerjee the singularities are of the following form

The singularity function can be given as:
q(x) = ,((⟨x-a⟩)^1)
An example of this is a beam  of unit length subjected by a load of the form shown in the figure
In this singularity function there is a slope discontinuity in q(x) at x = 1/3.

The second type of singularity function would be of the form:
q(x) = ((⟨x-a⟩)^2)
In this type of singularity function a beam of unit length is subjected to a load as shown in the accompanying figure. 
There is a discontinuity in the second derivative at x = 1/3.

These examples can be continued. In addition there are two very specific singularities:
	1) Concentrated force, and 2) Concentrated moment
1) A concentrated force is represented by 
	q(x) = ((⟨x-a⟩)^(-1))
In a more mathematical definition we can represent this force as a Dirac Delta function (δ)
	q(x) = δ(x-a)
This function is zero at any value other than a. At x = a, function becomes infinite and overall:
	(∫_(-∞)^∞)δ(x)dx = 1

2) It can be easily shown that a concentrated moment is represented by a derivative of a Delta function
	q(x) = ((⟨x-a⟩)^(-2))=δ'(x-a)

Solving Beam equation with mathematica subjected to various boundary conditions

We will now solve the beam equation subject to various boundary conditions. 
Consider a beam of length L = 1 subject to a concentrated load at x = a. The beam is fixed at both ends. This means that the displacement and the slope at x = 0 = 1 is zero. 

The beam equation for this problem will be given by:
E I  (((d^4)u)/(d (x^4))) = q(x) = δ(x-a)
Subject to boundary conditions:
u(0) = u(1) = 0;    u'(0) = u'(1) = 0. 

Without losing any generality we can say that EI = 1 and the length of the beam is one. We now solve this fourth order differential equation subject to these four boundary conditions using the function DSolve.

DSolve[{v''''[x]==DiracDelta[x-a], v[0]==0, v[1]==0, v'[1]==0, v'[0]==0}, v[x], x]//FullSimplify
u[x_, a_]= v[x]/.%[[1]];

The bending moment and the shear force are obtained by taking appropriate derivatives.

M[x_, a_]= D[u[x,a], x,x]//FullSimplify
V[x_, a_]= D[M[x,a], x]//FullSimplify

In order to understand what these loads signify we need to make a plot of the Bending Moment and Shear Force. We make it by the following for a = 1/3.

Plot[M[x,1/3], {x, 0, 1}, BaseStyle->{Medium, FontFamily->"Times"}, PlotStyle->{Red, Thick}, AxesLabel->{"x", "Moment(x)"}]
Plot[V[x,1/3], {x, 0, 1}, BaseStyle->{Medium, FontFamily->"Times"}, PlotStyle->{Red, Thick}, AxesLabel->{"x", "Shear(x)"}]

But this is unsatisfactory because this does not tell us as to at what location of the Concentrated force will we expect highest bending moment and shear force. In order to do that we use the following neat command called Manipulate  in Mathematica.

Moment[a_]:= Plot[M[x,a], {x, 0, 1}, BaseStyle->{Medium, FontFamily->"Times"}, PlotStyle->{Red, Thick}, AxesLabel->{"x", "Moment(x)"}, PlotRange->All]
Shear[a_]:= Plot[V[x,a], {x, 0, 1}, BaseStyle->{Medium, FontFamily->"Times"}, PlotStyle->{Red, Thick}, AxesLabel->{"x", "Shear(x)"}, PlotRange->All]
Manipulate[Moment[a], {{a,1/2}, 0.05, 0.95, 0.05}]
Manipulate[Shear[a], {{a, 1/2}, 0.05, 0.95, 0.05}]

We can change the boundary conditions and solve the problem.

How can we solve this same problem using the standard approach without using the singularity functions.

Clear["Global`*"]

The beam is divided into two parts. We have the following equations governing each section of the beam.

S1 = DSolve[{v1''''[x]==0,v1[0]==0, v1'[0]==0, v1'''[a]==V1, v1''[a]==M1},v1[x],x] 
S2 = DSolve[{v2''''[x]==0, v2'''[a]==V2, v2''[a]==M2, v2[1]==0, v2'[1]==0}, v2[x],x]
u1[x_,a_]= v1[x]/.S1[[1]];
u2[x_,a_]=v2[x]/.S2[[1]];

We now obtain the slope, displacement, and the shear force at the point of application of the force. Once we obtain that we use the continuity conditions at x = a to evaluate the four unknown constants

Slope1[x_,a_]= D[u1[x,a],x];
Slope2[x_,a_]= D[u2[x,a],x];
Mom1[x_,a_]=D[u1[x,a],x,x];
Mom2[x_,a_]= D[u2[x,a],x,x];
Shear1[x_,a_]= D[Mom1[x,a],x];
Shear2[x_,a_]= D[Mom2[x,a],x];

There are four unknown constants V1, V2, M1, and M2. We know that the moments, displacements, and slopes are continous across the point of application of the load in the beam. Also the sum of the shear forces is equal to the applied load.

eq[1] = V1-V2-1//FullSimplify
eq[2] = M2-M1//FullSimplify
eq[3] = Slope1[a,a]-Slope2[a,a]//FullSimplify
eq[4] = u1[a,a]-u2[a,a]//FullSimplify
Solution = Solve[{eq[1]==0, eq[2]==0, eq[3]==0, eq[4]==0}, {M1, M2, V1, V2}]//FullSimplify

We now form a piece-wise function using the above values.

u1[x,a]/.Solution[[1]]//FullSimplify;
u2[x,a]/.Solution[[1]]//FullSimplify;
Mom1[x,a]/.Solution[[1]]//FullSimplify
Mom2[x,a]/.Solution[[1]]//FullSimplify
Shear1[x, a]/.Solution[[1]]//FullSimplify
Shear2[x,a]/.Solution[[1]]//FullSimplify

u[x_,a_]:= (1/6) (((-1+a))^2) (x^2) (x+a (-3+2 x))/;0≤ x<a;
u[x_,a_]:= (1/6) (a^2) (((-1+x))^2) (a-3 x+2 a x)/; a≤ x≤ 1;
Moment[x_,a_]:=(((-1+a))^2) (x+a (-1+2 x))/;0≤ x<a;
Moment[x_,a_]:=(a^2) (2-3 x+a (-1+2 x))/; a≤ x≤ 1;
Shear[x_,a_]:= (((-1+a))^2) (1+2 a)/;0≤ x<a;
Shear[x_, a_]:= (a^2) (-3+2 a)/; a≤ x≤ 1;

Plot[Shear[x,1/3], {x, 0, 1}]