Last week I created #14826 to solve statically indeterminate beam systems using boundary conditions (bc_slope and bc_deflection). This week I mainly focussed on implementing a logic to find the maximum bending moment and maximum shear force in a beam.

Initially, I thought it would be as simple as differentiating the bending_moment and shear force and then solving those using solve. But solve couldn’t represent Interval solutions and hence gave a NonImplemented error, as both of these quantities can occur in Intervals.

So instead of using solve over whole spam length, I found out points of discontinuity in the respective equations using

```
for term in terms:
if isinstance(term, Mul):
term = term.args[-1]
```

where terms are all the Muls extracted from Add. and `term.args[-1]`

gives us the point of singularity.

Now between two singularity points, our function can be:

- A continuous and differentiable function (hence no Interval solution)
- or a constant value

for the first scenario, you just use solve over that interval and see values at the endpoint. The higher one of both gives you maxima in that interval. For the second, the constant value is indeed maximum by itself. Then compare maxims of all intervals and return location and its value.

```
>>> from sympy.physics.continuum_mechanics.beam import Beam
>>> from sympy import symbols
>>> E, I = symbols('E, I')
>>> l, P = symbols('l, P', positive=True)
>>> b = Beam(l, E, I)
>>> R1, R2 = symbols('R1, R2')
>>> b.apply_load(R1, 0, -1)
>>> b.apply_load(R2, l, -1)
>>> b.apply_load(P, 0, 0, end=l)
>>> b.solve_for_reaction_loads(R1, R2)
>>> b.max_bmoment()
(l/2, P*l**2/8)
```

**Next Week**

**Next Week**

- Beam class gives ValueError if passed value contains unit. So I would focus on fixing it.
- Read relevant theory for implementation of 3dBeam class.