I recommend you to start from the lowpass prototype, things will be a lot easier that way. I don't know what is the purpose of Cps
in your 2nd picture, or if it needs to be joined at the end (as @jonk mentions in the comments), but I will ignore it since your 1st and last pictures show the response of an ideal bandpass (otherwise you'd be having a zero in there, somewhere).
So, first, covert back the bandpass to its lowpass prototype by replacing the series LC with an L, and the parallel LC with a C. Then, since you're already using a software capable of symbolic calculation, you can use Kirchhoff on the whole circuit. This is how it looks like in wxMaxima but, the syntax should not be too unfamiliar:
eqs : [v1*(1/R_S + 1/s/L_1) = vi/R_S + v2/s/L_1,
v2*(1/s/L_1 +1/s/L_2 +s*C_1) = v1/s/L_1 + v3/s/L_2,
v3*(1/s/L_2 +1/s/L_3 +s*C_2) = v2/s/L_2 + vo/s/L_3,
vo*(1/s/L_3 +1/R_L) = v3/s/L_3];
solve(eqs, [v1, v2, v3, vo]);
You'll get this little monster (raw solution):
$$H(s)=\dfrac{1}{{C_1} {C_2} {L_1} {L_2} {L_3} {{s}^{5}}+\left( {C_1} {C_2} {L_2} {L_3}+{C_1} {C_2} {L_1} {L_2}\right) {{s}^{4}}+\left( {L_1} \left( {C_2} {L_3}+{C_1} {L_3}+{C_1} {L_2}\right) +{L_2} \left( {C_2} {L_3}+{C_1} {C_2}\right) \right) {{s}^{3}}+\left( {C_2} {L_3}+{C_1} {L_3}+\left( {C_2}+{C_1}\right) {L_2}+\left( {C_2}+{C_1}\right) {L_1}\right) {{s}^{2}}+\left( {L_3}+{L_2}+{L_1}+{C_2}+{C_1}\right) s+2} \tag{1}$$
Sticking to simplicity, the lowpass prototype transfer function for a 5th order Butterworth is:
$$\begin{align}
G(s)&=\dfrac{R_L}{R_L+R_S}\dfrac{1}{\left(s^2+2\sin\left(\dfrac{1}{10}\pi\right)s+1\right)\left(s^2+2\sin\left(\dfrac{3}{10}\pi\right)s+1\right)(s+1)} \\
&=\dfrac{0.5}{s^5+3.2361s^4+5.2361s^3+5.2361s^2+3.2361s+1} \tag{2}
\end{align}$$
Again, for simplicity, I've considered \$R_S=R_L=1\$ but, be careful (\$\color{red}{!!!}\$): having \$R_S\neq R_L\$ needs to be taken into account now. You can choose whichever one you see fit to be your reference, then make the other one relative to it (e.g. for \$R_S=50,\; R_L=75\$ make \$R_S=1,\; R_L=1.5\$).
You have several options to solve this: either use @Verbal Kint's method, or filter synthesis, or brute-force approach. Since you already have the software, I'll choose the 3rd (it's also easier to type)
NB: a 5th order will push the limits of the numeric solver. Higher orders might not even be able to be calculated without having to wait a LOT of time (if the solutions even come out). This is why I'd recommend the first two. Though, to be honest, where high orders are involved, I'd go with the 2nd since it will always find a solution, even if it will take a few steps to get there.
For this, normalize \$H(s)\$ to the zeroth power of \$s\$, such that the numerator and the last term are unity, both, and all the powers of \$s\$ will have terms which can be equated to the ones from \$G(s)\$. This will form a system of 5 equations with 5 unknowns, which can be fed to a solver. In wxMaxima, you get these solutions:
[[L_1=0.64546,L_2=1.9978,L_3=0.59284,C_1=1.6605,C_2=1.5756],
[L_1=0.63662-0.019807*%i,L_2=0.002248*%i+2.0,L_3=0.017559*%i+0.59944,C_1=1.6482-0.030148*%i,C_2=0.030148*%i+1.5879],
[L_1=0.019807*%i+0.63662,L_2=2.0-0.002248*%i,L_3=0.59944-0.017559*%i,C_1=0.030148*%i+1.6482,C_2=1.5879-0.030148*%i],
[L_1=0.61693-0.026175*%i,L_2=1.3418*10^-6*%i+2.0022,L_3=0.026174*%i+0.61692,C_1=1.618-0.042466*%i,C_2=0.042466*%i+1.618],
[L_1=0.026175*%i+0.61693,L_2=2.0022-1.3418*10^-6*%i,L_3=0.61692-0.026174*%i,C_1=0.042466*%i+1.618,C_2=1.618-0.042466*%i],
[L_1=0.59944-0.017574*%i,L_2=2.0-0.0022499*%i,L_3=0.019824*%i+0.63662,C_1=1.5879-0.030173*%i,C_2=0.030173*%i+1.6482],
[L_1=0.017574*%i+0.59944,L_2=0.0022499*%i+2.0,L_3=0.63662-0.019824*%i,C_1=0.030173*%i+1.5879,C_2=1.6482-0.030173*%i],
[L_1=0.59283,L_2=1.9978,L_3=0.64548,C_1=1.5756,C_2=1.6605],
[L_1=0.64542,L_2=1.9978,L_3=0.59288,C_1=1.6604,C_2=1.5756],
[L_1=0.63644-0.01959*%i,L_2=0.0022026*%i+2.0,L_3=0.017387*%i+0.59962,C_1=1.6479-0.029836*%i,C_2=0.029836*%i+1.5882],
[L_1=0.01959*%i+0.63644,L_2=2.0-0.0022026*%i,L_3=0.59962-0.017387*%i,C_1=0.029836*%i+1.6479,C_2=1.5882-0.029836*%i],
[L_1=0.61694-0.026152*%i,L_2=2.7596*10^-6*%i+2.0022,L_3=0.026149*%i+0.61691,C_1=1.6181-0.042427*%i,C_2=0.042427*%i+1.618],
[L_1=0.026152*%i+0.61694,L_2=2.0022-2.7596*10^-6*%i,L_3=0.61691-0.026149*%i,C_1=0.042427*%i+1.6181,C_2=1.618-0.042427*%i],
[L_1=0.59962-0.017418*%i,L_2=2.0-0.0022065*%i,L_3=0.019625*%i+0.63644,C_1=1.5882-0.029889*%i,C_2=0.029889*%i+1.6479],
[L_1=0.017418*%i+0.59962,L_2=0.0022065*%i+2.0,L_3=0.63644-0.019625*%i,C_1=0.029889*%i+1.5882,C_2=1.6479-0.029889*%i],
[L_1=0.59284,L_2=1.9978,L_3=0.64546,C_1=1.5756,C_2=1.6605]]
They look intimidating until you realize that only the first, the last, and the two middle rows have real solutions, which means you can choose any one of those four and all of them will give a 5th order Butterworth. Below is a SPICE confirmation, where the four outputs have been slightly offset to better see that they are all there (their phases, though, completely overlap, confirming identical results):

All that's left to do now is to convert it back to a bandpass with:
$$s\rightarrow\dfrac{s^2+\omega_c^2}{BW s}=\dfrac{s}{BW}+\dfrac{\omega_c^2}{BW s}=\dfrac{s}{BW}+\dfrac{1}{s\frac{BW}{\omega_c^2}} \tag{3}$$
What (3) really shows is that the \$s\$ operator, signifying a series inductor (\$sL\$), is replaced by a series combination of an inductor and a capacitor; the opposite for the shunt elements. If the I/O resistances have different values, scale all L by whichever of the two you kept as reference (in the example above, by \$R_S\$) and C by its reciprocal.
As an example, to convert it to a bandpass with \$\omega_c=1\;\text{krad},\;BW=2\;\text{krad},\;R_S=R_L=50\;\Omega\$, this is the result:

The readings are only different in the decimals, because:
- the values have been rounded (some to no decimals)
- the inductors have a default
Rser=1m
(to avoid noobs hooting themselves in the foot when placing a voltage source in parallel with an inductor)
- the number of points/decade is "only"
101
, which might affect some decimals (maybe not the first one)
And the values were taken as the last row from the solutions, then calculated as follows:
- for
Lb1
: 0.59284*50/2e3
- for
Cb1
: 1/0.59284/50*2e3/1e6
- for
Lb2
: 1/1.5756*50*2e3/1e6
- for
Cb2
: 1.5756/50/2e3
- ...