How to combine two plots in same figure with axis break?

I need a figure like this

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is a program for 1-D Photonic crystal 
% For Photonic band structure calculation
% This program uses Transfer Matrix Method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc 
clear
close all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define material properties and simulation parameters
P = [0 50 100 150 200];          %%% Hydrostatic pressure in MPa
n0 = 1;         % Refractive index of air
ns = 1.46;      % Refractive index of subtrate
n01 = 1.578;    % Refractive index of PS
n02 = 1.484;    % Refractive index of PMMA
% Calculate Epsilon for n01, n02
e01 = n01^2;    e02 = n02^2;
% Constants
p11_1 = 0.32;   p12_1 = 0.31;   p11_2 = 0.3;    p12_2 = 0.297;
v1 = 0.35;
E1 = 3.3e3; % Convert to scientific notation (3.3 x 10^3)
v2 = 0.37;
E2 = 3.0303e3; % Convert to scientific notation (3.0303 x 10^3)
% Formula for calculating nA
% Initialization of parameters
dA=780e-9;                          % Thickness of First Layer in meters
dB=830e-9;                          % Thickness of Second Layer in meters
c=3e8;                              % Velocity of Light (m/s)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Frequency Range
f=linspace(58,67,16001);                    % Frequncy in THz
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k = 1:length(P)
for loop=1:length(f)
    
e1 = e01 - ((e01^2)/2) * (((p11_1/E1) * (v1+1) * P(k)) + ((p12_1/E1) * (3*v1+1) * P(k)));
nA = sqrt(e1);
% Formula for calculating nA
e2 = e02 - ((e02^2)/2) * (((p11_2/E2) * (v2+1) * P(k)) + ((p12_2/E2) * (3*v2+1) * P(k)));
nB = sqrt(e2);
    
    %%% angular frequency
    
    w=2*pi*f*1e12;
    a=dA+dB;                                         %%% Lattice Constant   
   DA=((w(loop))/c)*dA*nA;                          %%% wave number of first layer
   DB=((w(loop))/c)*dB*nB;                          %%% wave number of second layer   
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   %%% Transfer Matrix elements of first layer
   m11=cos(DA); m12=-1i*sin(DA)/nA; m21=-1i*nA*sin(DA); m22=cos(DA);
   mA=[m11 m12; m21 m22];
   %%% Transfer MAtrix elements of Second layer
   l11=cos(DB); l12=-1i*sin(DB)/nB; l21=-1i*nB*sin(DB); l22=cos(DB);
   mB=[l11 l12; l21 l22];
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%% Unit cell   
   m=(mA*mB);
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%% Photonic Band Structure
   K(loop,k)=(1/(pi))*acos((m(1,1)+m(2,2))/2);
end 
   
end
figure(1)
subplot(1,2,1)
plot(-real(K),f)

legend('P=0','P=50','P=100','P=150','P=200')
% Set new axis limits to cut the x-axis between -0.9 and -0.5
xlim([-1 -0.9])
ylabel("Frequency (THz)")
xlabel("Wave vector")
subplot(1,2,2)
plot(real(K),f)
legend('P=0','P=50','P=100','P=150','P=200')
xlim([0.9 1])
xlabel("Wave vector")

Answer 1

[Voss]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is a program for 1-D Photonic crystal 
% For Photonic band structure calculation
% This program uses Transfer Matrix Method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc 
clear
close all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define material properties and simulation parameters
P = [0 50 100 150 200];          %%% Hydrostatic pressure in MPa
n0 = 1;         % Refractive index of air
ns = 1.46;      % Refractive index of subtrate
n01 = 1.578;    % Refractive index of PS
n02 = 1.484;    % Refractive index of PMMA
% Calculate Epsilon for n01, n02
e01 = n01^2;    e02 = n02^2;
% Constants
p11_1 = 0.32;   p12_1 = 0.31;   p11_2 = 0.3;    p12_2 = 0.297;
v1 = 0.35;
E1 = 3.3e3; % Convert to scientific notation (3.3 x 10^3)
v2 = 0.37;
E2 = 3.0303e3; % Convert to scientific notation (3.0303 x 10^3)
% Formula for calculating nA
% Initialization of parameters
dA=780e-9;                          % Thickness of First Layer in meters
dB=830e-9;                          % Thickness of Second Layer in meters
c=3e8;                              % Velocity of Light (m/s)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Frequency Range
f=linspace(58,67,16001);                    % Frequncy in THz
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k = 1:length(P)
    for loop=1:length(f)
        
        e1 = e01 - ((e01^2)/2) * (((p11_1/E1) * (v1+1) * P(k)) + ((p12_1/E1) * (3*v1+1) * P(k)));
        nA = sqrt(e1);
        % Formula for calculating nA
        e2 = e02 - ((e02^2)/2) * (((p11_2/E2) * (v2+1) * P(k)) + ((p12_2/E2) * (3*v2+1) * P(k)));
        nB = sqrt(e2);
        
        %%% angular frequency
        
        w=2*pi*f*1e12;
        a=dA+dB;                                         %%% Lattice Constant   
        DA=((w(loop))/c)*dA*nA;                          %%% wave number of first layer
        DB=((w(loop))/c)*dB*nB;                          %%% wave number of second layer   
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        %%% Transfer Matrix elements of first layer
        m11=cos(DA); m12=-1i*sin(DA)/nA; m21=-1i*nA*sin(DA); m22=cos(DA);
        mA=[m11 m12; m21 m22];
        %%% Transfer MAtrix elements of Second layer
        l11=cos(DB); l12=-1i*sin(DB)/nB; l21=-1i*nB*sin(DB); l22=cos(DB);
        mB=[l11 l12; l21 l22];
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%% Unit cell   
        m=(mA*mB);
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%% Photonic Band Structure
        K(loop,k)=(1/(pi))*acos((m(1,1)+m(2,2))/2);
    end 
end

figure()

x = real([-K K]); % data
xtl = [-1 -0.95 -0.9 0.9 0.95 1]; % x-tick labels
xl = [-1 -0.87; 0.87 1]; % x-limits
gap_size = 0.05;

% transform data to (0,1) and split into two parts based on xl
[xt1,xt2] = x_transform(x,xl,gap_size);

% use 5 colors:
c = colororder();
colororder(c(1:size(xt1,2),:))

% plot 10 lines:
plot(xt1,f)
hold on
plot(xt2,f)

ylabel("Frequency (THz)")
xlabel("Wave vector")

ax = gca();

% legend:
lgd = legend('P=0','P=50','P=100','P=150','P=200');
title(lgd,'Hydrostatic Pressure')
lgd.EdgeColor = 'none';
lgd.FontWeight = 'bold';
lgd.Position(1) = ax.Position(1)+(ax.Position(3)-lgd.Position(3))/2;
lgd.Position(2) = ax.Position(2)+ax.Position(4)/2-lgd.Position(4);

% axes appearance:
[xtlt1,xtlt2] = x_transform(xtl,xl,gap_size);
set(ax, ...
    'FontWeight','bold', ...
    'Box','on', ...
    'XLim',[0 1], ...
    'XTick',[xtlt1 xtlt2], ...
    'XTickLabel',xtl);

% gaps and gap lines:
yl = get(ax,'YLim');
dyl = yl(2)-yl(1);
set(ax,'YLimMode','manual');
patch( ...
    'Parent',ax, ...
    'XData',0.5+gap_size/2*[-1 -1; 0 0; 1 1; 0 0], ...
    'YData',yl+dyl*[1 -1].*[-1; 1; 1; -1]/30, ...
    'FaceColor','w', ...
    'EdgeColor','none', ...
    'Clipping','off', ...
    'HandleVisibility','off');
line( ...
    'Parent',ax, ...
    'XData',reshape(0.5+gap_size/2*[1 -1].*[-1; 0; NaN; 1; 0; NaN],1,[]), ...
    'YData',reshape(yl+dyl*[1 -1].*[-1; 1; NaN; 1; -1; NaN]/30,1,[]), ...
    'Color','k', ...
    'LineWidth',1, ...
    'Clipping','off', ...
    'HandleVisibility','off');

function [x1,x2] = x_transform(x,xl,gap_size)
[mx,nx] = size(x);
x1 = NaN(mx,nx);
x2 = NaN(mx,nx);
idx1 = x <= xl(1,2);
idx2 = x >= xl(2,1);
w = (1-gap_size)/2;
x1(idx1) = (x(idx1)-xl(1,1))/(xl(1,2)-xl(1,1))*w;
x2(idx2) = (x(idx2)-xl(2,1))/(xl(2,2)-xl(2,1))*w+1-w;
x1(:,all(isnan(x1),1)) = [];
x2(:,all(isnan(x2),1)) = [];
end