Draw decision boundary using Polynomial Kernel Function in SVM Matlab

I am trying to draw a decision boundary for fisherisis dataset on Matlab. Below is the code that I am working on:

load fisheriris;
X = meas(:,1:2);
y = species;

% PLOT
figure;
plot(X(strcmp(y,"setosa"),1),X(strcmp(y,"setosa"),2),'.','DisplayName','setosa'); hold on
plot(X(strcmp(y,"versicolor"),1),X(strcmp(y,"versicolor"),2),'.','DisplayName','versicolor');
plot(X(strcmp(y,"virginica"),1),X(strcmp(y,"virginica"),2),'.','DisplayName','virginica'); hold off
legend("setosa","versicolor","virginica")

% SVM
t = templateSVM("KernelFunction","polynomial", "BoxConstraint",1,"KernelScale","auto","PolynomialOrder",2,"Standardize",1);
SVMModel = fitcecoc(X,y,"Learners", t);

% DECISION BOUNDARY
%%% Draw decision boundary here

The above SVMModel I am using is classified multiclass SVM (3 classes: "setosa", "versicolor", "virginica") and 2 first column of the dataset for the train data.

I am expecting that I can plot 2 or 3 lines or curves displaying the boundary of each species in my SVMModel. I have read several websites that somewhat included my problems but nothing works for me because they are using linear SVM or they are using linear equation... Below is list of websites that I have searched about:

1. First website, but this is for binary classification, in my case, I have 3 classes.
2. Second website, but this is for linear SVM, in my case is polynomial kernel function.
3. Third website, but in this situation, the properties of ecoc has the property Beta, in my case the Beta is not provided.

Anyone have any idea how to draw decision boundary in my situation? Thanks a lot in advance.

Add these lines after your existing code:

[x1Grid,x2Grid] = meshgrid(linspace(min(X(:,1)), max(X(:,1)), 100), ...
                          linspace(min(X(:,2)), max(X(:,2)), 100));
[~,scores] = predict(SVMModel, [x1Grid(:), x2Grid(:)]);
scores = reshape(scores, [size(x1Grid,1), size(x1Grid,2), 3]);

% add color to each region
[pred,~] = predict(SVMModel, [x1Grid(:), x2Grid(:)]);
contourf(x1Grid, x2Grid, reshape(double(categorical(pred)), size(x1Grid)), 'LineStyle', 'none')
contourcmap('jet', 'Colorbar', 'on')

contour(x1Grid, x2Grid, scores(:,:,1) - max(scores(:,:,2:3),[],3), [0 0], 'k-');
contour(x1Grid, x2Grid, scores(:,:,2) - max(scores(:,:,[1 3]),[],3), [0 0], 'k--');
contour(x1Grid, x2Grid, scores(:,:,3) - max(scores(:,:,1:2),[],3), [0 0], 'k:');

The outcome is like: