Python: How to fit data to a model (lmfit) that includes multiple breakpoints (discrete values)?

Question

I would like to fit parameters p1-p6 into the function f(x1, x2) with multiple discrete conditions - the math would look like that:

if (x1 > p4) & (x2 < p6):
  y = p1*x1
elif x1 > p5:
  y = p2*x2
else:
  y = p3 * np.maximum(x1, x2)

I know how to solve the problem if I do not have this part of the equation: & (x2 < p6) so something like this:

if x1 > p4:
  y = p1*x1
elif x1 > p5:
  y = p2*x2
else:
  y = p3 * np.maximum(x1, x2)

in such case I would have to convert the discrete part into the continuous using scipy.special.erf():

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from lmfit import Model, minimize, report_fit, Parameters
from scipy.special import erf # to convert discrete values

def if_else_model(x, p1, p2, p3, p4, p5):
    
    # I have two arguments in x:
    x_1 = x[0]
    x_2 = x[1]
     
    # I have three functions:
    y1 = p1*x_1
    y2 = p2*x_2 * np.ones(len(x_2))
    y3 = p3*np.maximum(x_1, x_2) * np.ones(len(x_2))
    
    # I change discrete to continuous by adding errror: (erf(x-parameter)+1)/2
    # And as I have if - elif - else: I used y_2 to split it to: if - else - if - else - I guess there is a  better way?
    y_2 = y2 * (erf(x_1-p5)+1)/2 + y3 * (1-erf(x_1-p5))/2    
    y = y1 * (erf(x_1-p4)+1)/2 + y_2 * (1-erf(x_1-p4))/2
  
    return y

# I add the noise to create fake data to fit the model:
def if_else_model_fakedata(x, p1, p2, p3, p4, p5):
    return if_else_model(x, p1, p2, p3, p4, p5) + np.random.normal(size=x.size, scale=1.2)+4

# Define my x data:
x1 = np.linspace(1,50,50)
x2 = np.linspace(10,20,50)
x = [x1, x2]`

# Create the fake data
y = if_else_model_fakedata(x2,2, 2, 5, 30, 10)

# We can visualize it based on one argument (x1):
plt.plot(x1,y,"ko")
plt.plot(x1, if_else_model(x, 1, 0.2, 2, 25, 15), "r--")

# Model definition and initiating parameters:
test_model = Model(if_else_model)
params = test_model.make_params(p1=1, p2 = 0.2, p3 = 2, p4 = 25, p5 = 15)

# Fitting
result = test_model.fit(y, params, x=x)

# Plotting results
plt.plot(x1,y,"ko")
plt.plot(x1,result.init_fit, "r--")
plt.plot(x1,result.best_fit, "g-")

Thus, again the question:

how to add another condition for the second part--> & (x2 < p6)?

I am lost with converting discrete conditions into continuous numbers - step by step guidance would be great as I'm a Python beginner.