Multivariate time series forecast with VAR confusion
In essence, you want to both incorporate the past historical values of the target timeseries and the (past and) current historical values of other timeseries to predict the current value of the target timeseries.
In the tutorial you gave, they define the VAR[1] model with two time series ($Y_{1, t}$, $Y_{2, t}$) as:
$$ Y_{1,t} = \alpha_{1} + \beta_{11,1}Y_{1,t-1} + \beta_{12,1}Y_{2,t-1} +\epsilon_{1,t}$$
$$ Y_{2,t} = \alpha_{2} + \beta_{21,1}Y_{1,t-1} + \beta_{22,1}Y_{2,t-1} +\epsilon_{2,t}$$
where the author defined $\alpha$ as the intercept, $\beta$ the lag coefficients and $\epsilon$ as the error term.
Let us consider $Y_{1,t}$ as the target time series. It seems that it could be modified such that it incorporates the past and current values of $Y_{2,t}$ as follows:
$$ Y_{1,t} = \alpha_{1} + \beta_{11,1}Y_{1,t-1} + \beta_{12,1}Y_{2,t-1} + \beta_{12,0}Y_{2, t} + \epsilon_{1,t}$$
Here, the term $\beta_{12,0}Y_{2, t}$ has been added, which signifies the dependence on the current value of the other time series. In fact, this equation is the only relevant one when determining the $\beta$ parameters, as $Y_{2,t}$ is known in the test set.
Generalizing this to an arbitrary number of time series $n$ and lag $l$, the model for your target timeseries becomes:
$$ Y_{1,t} = \text{VAR}[l] + \sum^{n}_{i=2}\beta_{1i,0}Y_{i, t}$$
with VAR$[l]$ indicating the original VAR model from the article.
EDIT:
Worked example based on sample data you provided. We have time series $X,Y,Z$. Let us model our time series with lag = 1. Then our model for $Z_{t}$ becomes:
$$Z_{t} = (\alpha + \beta_{ZZ, 1}Z_{t-1} + \beta_{ZY,1}Y_{t-1} + \beta_{ZX,1}X_{t-1} + \epsilon_{t}) + \beta_{ZY,0}Y_{t} + \beta_{ZX,0}X_{t}$$
The values between brackets is standard vector autoregression. We can then use a simple linear regression procedure. Note that we only have two samples to use for fitting, because at Day 1 there are no previous timepoints. Thus we have the following equations:
$$Z_{2} = 60 = (\alpha + \beta_{ZZ, 1}Z_{1}+ \beta_{ZY,1}Y_{1} + \beta_{ZX,1}X_{1}) + \beta_{ZY,0}Y_{2} + \beta_{ZX,0}X_{2}) = (\alpha + \beta_{ZZ, 1}30 + \beta_{ZY,1}20 + \beta_{ZX,1}11) + \beta_{ZY,0}40 + \beta_{ZX,0}22)$$
$$Z_{3} = 90 = (\alpha + \beta_{ZZ, 1}Z_{2}+ \beta_{ZY,1}Y_{2} + \beta_{ZX,1}X_{2}) + \beta_{ZY,0}Y_{3} + \beta_{ZX,0}X_{3}) = (\alpha + \beta_{ZZ, 1}60 + \beta_{ZY,1}40 + \beta_{ZX,1}22) + \beta_{ZY,0}60 + \beta_{ZX,0}33)$$
Using python we can solve this as follows:
from sklearn.linear_model import LinearRegression

#Z-vals
Z_1 = 30
Z_2 = 60
Z_3 = 90

#Y-vals
Y_1 = 20
Y_2 = 40
Y_3 = 60

#X-vals
X_1 = 11
X_2 = 22
X_3 = 33

LR = LinearRegression()

result = LR.fit([[Z_1, Y_1, X_1, Y_2, X_2],
                [Z_2, Y_2, X_2, Y_3, X_3]],[Z_2, Z_3])

print(f'Alfa is {result.intercept_}')
print(f'Parameters are {result.coef_}')


which yields
Alfa is 21.95159629248201
Parameters are [0.46343975 0.30895984 0.16992791 0.30895984 0.16992791]

And thus our predictive model for Z has become (after rounding):
$$Z_{t} = 21.95 + 0.46Z_{t-1} + 0.31Y_{t-1} + 0.17X_{t-1} + 0.31Y_{t} + 0.17X_{t}$$